text:
# The antiferromagnetic Ising chain in a mixed transverse and longitudinal magnetic field

## I Introduction

Recently, the study of the field-induced effects in low-dimensional quantum spin systems has been attracting much interest from theoretical and experimental points of view YbAs ,Affleck ,kufo ,Kohgi . For the system with an anisotropy of exchange interactions the magnetic properties essentially depend on the direction of the applied magnetic field. For example, the behavior of the one-dimensional antiferromagnetic $`XXZ`$ model in a transverse magnetic field is drastically different in comparison with the case of the longitudinal field. In particular, the transverse field induces the staggered magnetization in the perpendicular direction and the continuous phase transition takes places at some critical field mori ,hieda ,we ,dutta . This effect has been observed in quasi-one-dimensional antiferromagnet $`\mathrm{Cs}_2\mathrm{CoCl}_4`$ kenz , where the magnetic field has both the transverse and longitudinal components. Therefore, it is important to study the properties of the antiferromagnetic $`s=\frac{1}{2}`$ $`XXZ`$ model in mixed transverse and longitudinal magnetic fields

$$H=(S_n^xS_{n+1}^x+S_n^yS_{n+1}^y+\mathrm{\Delta }S_n^zS_{n+1}^z)h_xS_n^xh_zS_n^z$$
(1)
We consider the most simple case of this model – the antiferromagnetic Ising chain given by the Hamiltonian

$$H=S_n^zS_{n+1}^zh_xS_n^xh_zS_n^z$$
(2)
In spite of the simple form of this Hamiltonian it cannot be solved exactly. In sen  the ground state phase diagram of this model has been investigated using numerical diagonalization of finite systems and finite size scaling procedure. It was found sen  that the transition line between the ordered and disordered phases exists and it was assumed that the model belongs to the universality class of the two-dimensional Ising model.

We note that the ferromagnetic Ising chain in the mixed fields has been studied intensively fog ,zam ,fateev . Though the ferromagnetic and the antiferromagnetic Ising chains in the transverse magnetic field are equivalent, properties of these two models are very different at $`h_z0`$. (In fact, the model (2) can be transformed to the ferromagnetic chain but in a staggered longitudinal field.) For example, the ground state phase transition in the ferromagnetic model is smeared out by the longitudinal field in contrast to the antiferromagnetic model for which the phase transition remains at $`h_z0`$.

In this paper we study the model (2) using the Density Matrix Renormalization Group (DMRG) technique white . This method allows us to consider the systems up to a few hundred sites and to determine the transition line with high accuracy. Using the finite-size estimation of the ground state energy and low-lying excitations we will show that the model (2) on the transition line is described by the conformal field theory with the central charge $`c=\frac{1}{2}`$. Besides, we consider the behavior of the model in the vicinity of the special points $`h_x=\frac{1}{2}`$, $`h_z=0`$ and $`h_x=0`$, $`h_z=1`$ and determine the form of the transition line near these points.

The paper is organized as follows. In Sec.II we will provide a qualitative physical picture of the ground state phase diagram based on the classical approach. In Sec.III and Sec.IV the behavior of the system in the vicinity of the special points of the critical line will be considered. In Sec.V we will present the DMRG calculation of the critical line. Sec.VI is devoted to the connection of the model (2) with the statistical quasi-one-dimensional Ising model in the external magnetic field. In Summary we discuss our results.

## II The classical approach

In order to provide a physical picture of the phase diagram of the model (2) we use the classical approximation, when spins are represented as three-dimensional vectors. The classical ground state is given by a configuration in which all spin vectors lie in the $`XZ`$ plane with the spins on odd and even sites pointing respectively at angles $`\phi _1`$ and $`\phi _2`$ with respect to the $`X`$ axis. The classical energy of this state is

$$E/N=\frac{1}{4}\mathrm{sin}\phi _1\mathrm{sin}\phi _2\frac{h_x}{4}(\mathrm{cos}\phi _1+\mathrm{cos}\phi _2)\frac{h_z}{4}(\mathrm{sin}\phi _1\mathrm{sin}\phi _2)$$
(3)
The angles $`\phi _1`$ and $`\phi _2`$, minimizing this energy are solutions of the following equations

$`\mathrm{cos}\phi _1\mathrm{sin}\phi _2+h_z\mathrm{cos}\phi _1h_x\mathrm{sin}\phi _1`$ $`=`$ $`0`$
$`\mathrm{sin}\phi _1\mathrm{cos}\phi _2h_z\mathrm{cos}\phi _2h_x\mathrm{sin}\phi _2`$ $`=`$ $`0`$ (4)
The solution of these equations is simple in the particular cases $`h_z=0`$

$`\phi _1`$ $`=`$ $`\phi _2=\mathrm{arccos}(h_x),h_x1`$
$`\phi _1`$ $`=`$ $`\phi _2=0,h_x>1`$ (5)
and $`h_x=0`$

$`\phi _1`$ $`=`$ $`\phi _2={\displaystyle \frac{\pi }{2}},h_z1`$
$`\phi _1`$ $`=`$ $`\phi _2={\displaystyle \frac{\pi }{2}},h_z>1`$ (6)
But in general case, when $`h_x0`$ and $`h_z0`$, the phase diagram is divided on two regions (see Fig. 1). In the paramagnetic region (PM) the energy minimum is given by the configuration with $`\phi _1=\phi _2=\phi `$, with $`\phi `$ determined by the equation

$$h_z\mathrm{cos}\phi h_x\mathrm{sin}\phi =\mathrm{sin}\phi \mathrm{cos}\phi $$
(7)
The antiferromagnetic (AF) region is characterized by non-zero staggered magnetizations in both $`X`$ and $`Z`$ directions

$`S_{2n}^xS_{2n+1}^x`$ $`=`$ $`\mathrm{sin}({\displaystyle \frac{\phi _2\phi _1}{2}})\mathrm{sin}({\displaystyle \frac{\phi _1+\phi _2}{2}})`$
$`S_{2n}^zS_{2n+1}^z`$ $`=`$ $`\mathrm{cos}({\displaystyle \frac{\phi _2\phi _1}{2}})\mathrm{sin}({\displaystyle \frac{\phi _1+\phi _2}{2}})`$ (8)
with $`\phi _1+\phi _20`$.

Thus, in the AF region the magnetic fields induce a perpendicular antiferromagnetic long-range order (LRO). In the classical approach the value $`(\phi _1+\phi _2)`$ plays the role of the LRO and it vanishes on the transition line determined by the equations

$`h_x\mathrm{cos}\phi +h_z\mathrm{sin}\phi `$ $`=`$ $`1`$
$`h_z\mathrm{cos}\phi h_x\mathrm{sin}\phi `$ $`=`$ $`\mathrm{sin}\phi \mathrm{cos}\phi `$ (9)
where $`\phi _1=\phi _2=\phi .`$

The solution of Eqs.(9) gives the transition line in the explicit form

$$h_z=\sqrt{1h_x^{2/3}}(1+h_x^{2/3})$$
(10)
The classical phase diagram is shown on Fig. 1. The transition (critical) line separates the phase with non-zero AFM LRO from the phase with uniform magnetization (the paramagnetic phase).

Of course, the classical approach does not give the correct description of the phase transition. At first, the quantum fluctuations shift the critical point $`h_x=1`$ to $`h_x=\frac{1}{2}`$ at $`h_z=0`$. Secondly, the form of the critical line is certainly incorrect at $`h_x1`$. Besides, the order parameter $`(\phi _1+\phi _2)`$ vanishes on the critical line with the critical exponent $`\frac{1}{2}`$, which is not valid at least in the critical point $`h_x=\frac{1}{2},h_z=0`$. Nevertheless, the fact of the generation of the staggered magnetizations perpendicular to the field direction is qualitatively true.

## III The critical ITF point

Before the numerical determination of the critical line we consider the behavior of the model (2) in the vicinity of two special points: $`h_x=\frac{1}{2},h_z=0`$ and $`h_x=0,h_z=1`$. This will give us the form of the critical line near these points.

In the case $`h_z=0`$ the model (2) reduces to the exactly solvable one-dimensional Ising model in the transverse field (the ITF model). This model is well studied Pfeuty  and it belongs to the universality class of the two-dimensional Ising model. At $`h_x<\frac{1}{2}`$ the ITF model is gapped and there is non-zero staggered magnetization $`(1)^nS_n^z`$. The ground state at $`h_x<\frac{1}{2}`$ is two-fold degenerated in the thermodynamic limit. At the point $`h_x=\frac{1}{2}`$ the model becomes gapless and the staggered magnetization vanishes with the critical exponent $`\frac{1}{8}`$. This is the transition point from the antiferromagnetic state to the paramagnetic gapped state. The gap vanishes linearly at $`h_x=\frac{1}{2}`$.

Now we consider the ITF model with small longitudinal magnetic field. For small magnetic field $`h_z`$ we rewrite the Hamiltonian (2) in the form

$`H`$ $`=`$ $`H_0+V`$
$`H_0`$ $`=`$ $`{\displaystyle S_n^zS_{n+1}^z}h_x{\displaystyle S_n^x}`$
$`V`$ $`=`$ $`h_z{\displaystyle S_n^z}`$ (11)
After Jordan-Wigner transformation to the Fermi operators $`c_n^+`$, $`c_n`$

$`S_n^z`$ $`=`$ $`{\displaystyle \frac{c_n^++c_n}{2}}{\displaystyle \underset{j<n}{}}\left(12c_j^+c_j\right)`$
$`S_n^x`$ $`=`$ $`{\displaystyle \frac{1}{2}}c_n^+c_n`$ (12)
the ITF Hamiltonian $`H_0`$ takes a bilinear form

$$H_0=h_x\left(c_n^+c_n\frac{1}{2}\right)+\frac{1}{4}\left(c_n^+c_n\right)\left(c_{n+1}^++c_{n+1}\right)$$
(13)
The Hamiltonian $`H_0`$ commutes with the ’parity’ operator

$$P=\mathrm{exp}\left(i\pi c_n^+c_n\right)$$
(14)
because $`H_0`$ can change the number of $`c_n`$ excitation by an even number only. Therefore, the space of states of $`H_0`$ is divided on two sectors with odd ($`P=1`$) and even ($`P=1`$) number of the Fermi particles $`c_n`$.

The Hamiltonian $`H_0`$ is diagonalized exactly LSM

$$H_0=\epsilon _k\left(\eta _k^+\eta _k\frac{1}{2}\right)$$
(15)
with Fermi particles $`\eta _k`$ and spectrum

$$\epsilon _k^2=h_x^2+\frac{1}{4}+h_x\mathrm{cos}k$$
(16)
which gives a gap at momentum $`\pi `$

$$m=\left|h_x\frac{1}{2}\right|$$
(17)
The gap vanishes at the critical point $`h_x=1/2`$, where the spectrum becomes

$$\epsilon _k=\mathrm{cos}\frac{k}{2}$$
(18)
Below in this section we will consider the perturbation theory in $`V`$ (11) in the critical point $`h_x=1/2`$.

The transition operator $`S^z=S_n^z`$ in $`V`$ (11) conserves the momentum and changes the ’parity’, because as it follows from Eq.(12) it changes the number of the Fermi particles $`c_i`$ (and also the number of $`\eta _k`$ particles) by odd number. Therefore, the non-zero matrix elements with the transition operator $`S^z`$ have the states with equal momentum but different parity. The last fact means that the perturbation theory in $`V`$ contains only even orders.

The second order correction to the ground state energy is:

$$\delta E_0^{(2)}=h_z^2\underset{s}{}\frac{0\left|S^z\right|ss\left|S^z\right|0}{E_0E_s}$$
(19)
The ground state $`|0`$ has the momentum $`q=0`$ and zero number of the $`\eta _k`$ particles ($`P=1`$). Therefore, the non-zero contribution to the sum in (19) is given by the intermediate states $`|s`$ with zero momentum and the parity $`P=1`$. As follows from Eq.(18), all states $`|s`$ with momentum $`q=0`$ and odd number of the $`\eta _k`$ particles ($`P=1`$) have ‘high’ energies $`E_sE_0\epsilon _s1`$. On the contrary, among the states with $`q=0`$ and even number of the $`\eta _k`$ particles ($`P=1`$) there are many states like $`\eta _{\pi k}^+\eta _{\pi +k}^+|0`$ with small $`k`$ having small excitation energies $`\epsilon _sk`$ and they can lead to infrared divergencies.

Hereinafter we consider large but finite systems of length $`N`$. We shall study the dependence of dominant contributions to the perturbation theory on $`N`$, omitting numerical factors. Using the fact that $`\epsilon _s1`$ one can rewrite Eq.(19) as

$$\delta E_0^{(2)}h_z^2\underset{s}{}0\left|S^z\right|ss\left|S^z\right|0=h_z^20\left|\left(S^z\right)^2\right|0$$
(20)
In higher orders of the perturbation series for the ground state energy each second intermediate state $`|s`$ has odd number of the Fermi particles ($`P=1`$), and, therefore, high energy $`\epsilon _s1`$. For example, let us consider the fourth order correction to the ground state energy

$`\delta E_0^{(4)}`$ $`=`$ $`h_z^4{\displaystyle \underset{s,s^{},s^{\prime \prime }}{}}{\displaystyle \frac{0\left|S^z\right|ss\left|S^z\right|s^{}s^{}\left|S^z\right|s^{\prime \prime }s^{\prime \prime }\left|S^z\right|0}{\left(E_0E_s\right)\left(E_0E_s^{}\right)\left(E_0E_{s^{\prime \prime }}\right)}}`$ (21)
$`\delta E_0^{(2)}h_z^2{\displaystyle \underset{s}{}}{\displaystyle \frac{0\left|S^z\right|ss\left|S^z\right|0}{\left(E_0E_s\right)^2}}`$
All intermediate states $`s,s^{},s^{\prime \prime }`$ have momentum $`q=0`$. The states $`s^{}`$ have $`P=1`$ and some of them have small excitation energies $`\epsilon _s^{}1/N`$, while the states $`s`$ and $`s^{\prime \prime }`$ have $`P=1`$ and high energies $`\epsilon _s1`$. Therefore, one can sum over the intermediate states $`s`$ and $`s^{\prime \prime }`$, which reduces Eq.(21) to

$$\delta E_0^{(4)}h_z^4\underset{s^{}}{}\frac{0\left|\left(S^z\right)^2\right|s^{}s^{}\left|\left(S^z\right)^2\right|0}{E_0E_s^{}}$$
(22)
Now we note, that this expression looks like the second order correction and Eq.(20) looks like the first order correction to the ground state energy with a perturbation $`h_z^2\left(S^z\right)^2`$. In a similar way one can sum over all intermediate states with $`q=0`$ and $`P=1`$ in all orders of perturbation series. As a result, we arrive at the perturbation theory with the effective perturbation

$$V_1=h_z^2\left(S^z\right)^2=h_z^2\underset{n,m}{}S_n^zS_m^z$$
(23)
We note, that the perturbation theory with the perturbation $`V_1`$ coincides with the original perturbation theory (11) in a sense, that both perturbation series have the same order of divergencies (or power of $`N`$) at each order in $`h_z`$. But numerical factors at each order in $`h_z`$ can be different.

The perturbation $`V_1`$ commutes with the parity operator (14) and conserves the momentum. Therefore, the perturbation series in $`V_1`$ contains the intermediate states with $`q=0`$ and $`P=1`$ only, and some of these states have small excitation energies $`\epsilon _s1/N`$. These states give the most divergent contribution to the perturbation series in $`V_1`$ and further we shall take into account these states only.

Now we need to estimate the matrix elements of the operator $`\left(S^z\right)^2`$. The behavior of the correlation function $`S_n^zS_m^z`$ in the ground state and in the low-lying states with excitation energies $`\epsilon _s1/N`$ on large distances is known McCoy

$$S_n^zS_m^z\frac{(1)^{nm}}{\left|nm\right|^{1/4}}$$
(24)
and, therefore, due to oscillation of the correlator $`S_n^zS_m^z`$ the sum over $`n`$ and $`m`$ can be estimated as

$$\underset{n,m}{}S_n^zS_m^z=\frac{N}{4}+2N\underset{n>1}{}S_1^zS_n^z0.07465(1)N$$
(25)
where the constant was found from extrapolation of the exact results for finite chains with $`N=6,\mathrm{}14`$.

The non-diagonal matrix elements of the operator $`\left(S^z\right)^2`$ with two different low-lying states was also calculated numerically. It was found that the only non-zero matrix elements (in the thermodynamic limit) are given by the states $`s,s^{}`$ differing by two $`\eta _k`$ particles, like $`|s^{}=\eta _k^+\eta _k^+|s`$. All such pairs of states give the same value for the matrix element

$$s\left|\left(S^z\right)^2\right|s^{}=0.3108(1)$$
(26)
while all other matrix elements exponentially drop with $`N`$.

Now we estimate all terms of perturbation series in $`h_z`$. According to Eq.(25) the second order correction (20) is proportional to $`h_z^2N`$. The factor at $`h_z^2N`$ was found numerically by exact diagonalization of finite systems and the following calculation of the sum (19). The data for $`N=6,\mathrm{}14`$ are well extrapolated and give

$$\delta E_0^{(2)}=0.07060(5)h_z^2N$$
(27)
According to Eq.(27) the zero-field susceptibility is $`\chi _z=0.1412(1)`$. This result is in a perfect agreement with the value $`\chi _z=0.14118\mathrm{}`$ obtained analytically in muller .

Using Eq.(26) we find that the non-zero contribution to the fourth order correction to the ground state energy (22) is given by the states $`|s=\eta _k^+\eta _k^+|0`$ with the excitation energy $`\epsilon _s=2\mathrm{cos}\frac{k}{2}`$. Thus, because of small excitation energies at $`k\pi `$ in denominator in (22) the fourth order correction to the ground state energy turns out to be proportional to $`N`$ (we omit here logarithmic corrections)

$$\delta E_0^{(4)}h_z^4\underset{s}{}\frac{1}{\epsilon _s}h_z^4N$$
(28)
The $`m`$-th order in $`V_1`$ is proportional to $`h_z^{2m}`$. The denominator of the $`m`$-th order contains $`(m1)`$ small excitation energies $`1/N`$ and all matrix elements in the numerator are of the order of unity. Therefore, the $`m`$-th order in $`V_1`$ diverges as $`h_z^{2m}N^{m1}`$. But this is not valid for odd orders in $`V_1`$. The analysis shows that all odd orders in $`V_1`$ diverge as $`h_z^{2m}N^{m2}`$. For example, the third order correction in $`V_1`$ (the sixth order in $`h_z`$) for the ground state energy does not diverge:

$$\delta E_0^{(6)}h_z^6N$$
(29)
So, the odd orders in $`V_1`$ give the next order corrections to the ground state energy and we omit it below.

Summarize all above, we arrive at the perturbation series in the form

$$\delta E_0=a_0h_z^2Nh_z^4N\underset{n=0}{\overset{\mathrm{}}{}}b_n\left(h_z^2N\right)^n$$
(30)
with $`a_0=0.07059\mathrm{}`$ (Eq.(27)) and unknown constants $`b_n`$. One can see that the first divergence appears only in the eight-th order in $`h_z`$, and it is very difficult to observe it numerically. The exact numerical calculations up to $`N=14`$ of the second, the fourth and the sixth order corrections in $`h_z`$ confirm the form of the series (30).

The sum in Eq.(30) forms the scaling function $`f_0\left(x\right)`$ of the scaling parameter $`x=h_z^2N`$. Thus, the ground state energy takes the form

$$\delta E_0=a_0h_z^2Nh_z^4Nf_0\left(x\right)$$
(31)
Since the ground state energy is proportional to $`N`$, the scaling function $`f_0\left(x\right)`$ has finite thermodynamic limit at $`x\mathrm{}`$. Thus, the leading term of the perturbation theory for the ground state energy is given by the second order and the divergent part of the perturbation theory gives the correction $`h_z^4`$.

The perturbation series for low-lying states has the same form as in Eq.(31), but each low-lying state has its own scaling function. Therefore, for the first excited state (with momentum $`q=\pi `$) one has

$$\delta E_\pi =a_\pi h_z^2Nh_z^4Nf_\pi \left(x\right)$$
(32)
So, the mass gap $`m=\delta E_\pi \delta E_0`$ appears as:

$$m=a_mh_z^2+g\left(x\right)h_z^2$$
(33)
where $`a_m=(a_0a_\pi )N`$ and the scaling function $`g\left(x\right)=x\left(f_0\left(x\right)f_\pi \left(x\right)\right)`$. Since the gap is finite, the scaling function $`g\left(x\right)`$ in the thermodynamic limit at $`x\mathrm{}`$ must tend to some finite limit

$$m=a_mh_z^2+g\left(\mathrm{}\right)h_z^2$$
(34)
From the last equation we see that the gap is proportional to $`h_z^2`$, but the factor at $`h_z^2`$ is given not only by the second order correction $`a_m`$ but by all collected divergent orders of the perturbation series. The numerical estimation of the second order correction to the gap gives

$$a_m=0.1875(2)$$
(35)
In order to find the factor at $`h_z^2`$ in mass gap we performed DMRG calculations (for details see Sec.V) of the model (11) for $`h_z=0,\mathrm{}0.3`$ and $`N=20,\mathrm{}300`$. The dependence of the gap on the scaling parameter $`x=h_z^2N`$ is shown on Fig. 2. One can see that the points with different $`N`$ and $`h_z`$ lie perfectly on one curve. It definitely manifests that the scaling parameter is $`h_z^2N`$. The DMRG calculation data give the mass gap

$$m=0.37(1)h_z^2$$
(36)
From this equation follows that the second order correction to the mass gap (35) gives approximately one half of a total gap, and another half is collected by all other divergent orders contained in the scaling function $`g\left(x\right)`$.

We note, that in contrast to Eq.(36) for the ferromagnetic ITF model the mass gap at the critical point $`h_x=\frac{1}{2}`$ is proportional to $`h_z^{8/15}`$ zam ,fateev .

### III.1 Mapping to the effective model

Let us return to the estimation of the diagonal element of the operator $`V_1`$ (25). Using the fact of the oscillation of the correlator (24) in the ground state and in the low-lying states we can rewrite the sum in (25) approximately as a one half of the first term:

$$\left(S^z\right)^2=\frac{N}{4}+2\underset{n<m}{}S_n^zS_m^z\frac{N}{4}+\underset{n}{}S_n^zS_{n+1}^z$$
(37)
The last equation suggests that the perturbation $`V_1`$ can be reduced to the operator $`V_2`$

$$V_1V_2=ah_z^2Nbh_z^2S_n^zS_{n+1}^z$$
(38)
with some constants $`a,b`$, which we will define later. In order to verify this assumption we compared the matrix elements of the operator $`V_2`$ with those of $`V_1`$. We have found that the dependence of the matrix elements of the operator $`V_2`$ on $`N`$ is in a full accord with Eqs.(25) and (26). That is, the diagonal matrix elements of $`V_2`$ are proportional to $`N`$ and the only non-zero (non-diagonal) matrix elements have the pairs of states differing by two $`\eta _k`$ particles

$$s\left|S_n^zS_{n+1}^z\right|s^{}=\frac{1}{2}$$
(39)
Thus, we arrive at the effective Hamiltonian

$$H_{\mathrm{eff}}=ah_z^2N+(1bh_z^2)S_n^zS_{n+1}^z\frac{1}{2}S_n^x$$
(40)
Again, the original model (11) is equivalent to the effective model (40) in a sense, that the perturbation series for both models have the same order of divergencies (or power of $`N`$) at each order of $`h_z`$. But numerical factors at each order of $`h_z`$ can be different.

The advantage of this mapping lies in the fact that the effective model (40) is of the ITF type and, therefore, is exactly solvable one. The spectrum of $`H_{\mathrm{eff}}`$ is

$$\epsilon _k^2=\frac{1+J^2}{4}+\frac{J}{2}\mathrm{cos}k$$
(41)
with $`J=1bh_z^2`$. So, according to the spectrum of the effective model (41) the magnetic field $`h_z`$ produces the gap in the spectrum

$$m=\frac{b}{2}h_z^2$$
(42)
As was noted above, apart from numerical factors the effective model (40) has the same form of perturbation series in $`h_z`$ as the original model (11). Therefore, we can check the form of the perturbation theory (30) by studying the perturbation series in $`h_z`$ of the exactly solvable model (40).

For large but finite systems the ground state energy (with $`q=0`$) and the first excited state (with $`q=\pi `$) of the Hamiltonian (40) are known exactly hamer :

$`E_0`$ $`=`$ $`ah_z^2N{\displaystyle \frac{1}{4}}{\displaystyle \underset{j=1}{\overset{N}{}}}\sqrt{g^2+4(1g)\mathrm{sin}^2\left({\displaystyle \frac{\pi j}{N}}{\displaystyle \frac{\pi }{2N}}\right)}`$
$`E_\pi `$ $`=`$ $`ah_z^2N+{\displaystyle \frac{g}{4}}{\displaystyle \frac{1}{4}}{\displaystyle \underset{j=1}{\overset{N1}{}}}\sqrt{g^2+4(1g)\mathrm{sin}^2\left({\displaystyle \frac{\pi j}{N}}\right)}`$ (43)
where $`g=bh_z^2`$.

The perturbation series for the ground state energy and the gap can be easily found by formal expansion of these expressions in small parameter $`g`$. This results in the scaling form of perturbation series for the mass gap:

$$m=\frac{g}{4}+g^2N\underset{n=0}{\overset{\mathrm{}}{}}d_n\left(gN\right)^{2n}$$
(44)
where constants $`d_n`$ determine the scaling function $`G(x)`$ of a scaling parameter $`x=gN`$

$$m=\frac{g}{4}+gG(x)$$
(45)
The scaling function $`G(x)`$ was found in hamer . In the thermodynamic limit $`x\mathrm{}`$ the function $`G(x)1/4`$, resulting in the gap (42). One can see that the scaling parameter and the form of the perturbation series (44) coincides with Eq.(33).

Now we define the constants $`a`$ and $`b`$ in Eq.(40), so that the second order corrections to the ground state energy and the gap of the effective model (40) numerically coincide with Eqs.(27) and (35). This leads to the values

$`a`$ $`=`$ $`0.190(1)`$
$`b`$ $`=`$ $`0.750(1)`$ (46)
The mass gap of the effective model (40) with constants $`a`$ and $`b`$ defined above is

$$m=0.375(1)h_z^2$$
(47)
Surprisingly, the value of the mass gap of the effective model turns out to be very close to one found numerically for the original model (Eq.(36)). It means, that the gap scaling function $`G(x)`$ of the effective model (40) has the same thermodynamic limit as that for the model (11). Moreover, we have calculated numerically the fourth order corrections in $`h_z`$ to the ground state energy and the gap for the models (11) and (40) and found perfect numerical agreement (the second order corrections in $`h_z`$ for both models coincide by the definition (46)). Therefore, we expect that the mapping of the model (11) at $`h_x=1/2`$ to the model (40) with constants $`a`$ and $`b`$ defined in (46) is exact for low-lying excitations in the thermodynamic limit.

The mapping of the model (11) to the effective ITF model can be extended to the case $`h_x1/2`$. For this case the effective Hamiltonian becomes

$$H_{\mathrm{eff}}=a(h_x)h_z^2N+\left(1b(h_x)h_z^2\right)S_n^zS_{n+1}^zh_xS_n^x$$
(48)
where now the constants $`a,b`$ are the functions of $`h_x`$. The functions $`a(h_x)`$ and $`b(h_x)`$ are defined in such a way that the second order corrections to the ground state energy and the gap for the effective model (40) coincide with those for the original model (11).

Using the effective Hamiltonian (48) one can calculate the susceptibility $`\chi _z(h_x)`$. In the vicinity of the critical point $`h_x=1/2`$ it has a logarithmic singularity

$$\chi _z(h_x)=\chi _{zc}\frac{b}{\pi }\left(h_x\frac{1}{2}\right)\mathrm{log}\left(h_x\frac{1}{2}\right)$$
(49)
with $`\chi _{zc}=0.1412`$ from Eq.(27) and $`b=0.75`$ (Eq.(46)).

## IV The multicritical point

Another exactly solvable limit of the model (2) is the case $`h_x=0`$. The model at $`h_x=0`$ is the classical one. At $`h_z<1`$ its ground state is the antiferromagnet. At $`h_z=1`$ the first order phase transition to the ferromagnetic ground state occurs. This is so-called multicritical point. The Hamiltonian (2) at $`h_z=1`$ has a form

$`H`$ $`=`$ $`H_0+V`$
$`H_0`$ $`=`$ $`{\displaystyle \frac{N}{4}}+{\displaystyle \left(S_n^z\frac{1}{2}\right)\left(S_{n+1}^z\frac{1}{2}\right)}`$
$`V`$ $`=`$ $`{\displaystyle \frac{h_x}{2}}\left(S^++S^{}\right)`$ (50)
where $`S^\pm =S_n^\pm `$.

The ground state of $`H_0`$ is macroscopic degenerate: all spin configurations, excluding those with two neighbor spins pointing down, have the same energy $`\frac{N}{4}`$. The number of these states is $`(\frac{1+\sqrt{5}}{2})^N`$ domb . The transverse field $`h_x`$ lifts the degeneracy. The exact calculation of the first order correction in $`h_x`$ for $`N1`$ is rather complicated because it involves the exponentially large number of degenerate states. We carried out the approximate calculation of the perturbation theory in $`V`$ within degenerate manifold using a simple variational function in the form

$$\mathrm{\Psi }=\underset{m=0}{\overset{N/2}{}}c_m\mathrm{\Psi }_m$$
(51)
where $`\mathrm{\Psi }_m`$ are the sum (with equal weights) of all admissible states with $`m`$ spins down:

$$\mathrm{\Psi }_m=w_m^{\frac{1}{2}}(S^{}P)^m\mathrm{\Psi }_F$$
(52)
Here $`P`$ is a projector excluding states with two neighbor down spins, $`\mathrm{\Psi }_F`$ is the ferromagnetic state with all spins up and the normalization factors $`w_m`$ are

$$w_m=\frac{(Nm1)!N}{m!(N2m)!}$$
(53)
Matrix elements of $`V`$ with respect to $`\mathrm{\Psi }_m`$ are

$$\mathrm{\Psi }_mV\mathrm{\Psi }_m^{}=\frac{h_x}{2}(s_m\delta _{m^{},m1}+s_{m+1}\delta _{m^{},m+1})$$
(54)
where

$$s_m=\sqrt{\frac{m(N2m+2)(N2m+1)}{Nm}}$$
(55)
Coefficients $`c_m`$ in (51) obey the equations

$$(E+\frac{N}{4})c_m=\frac{h_x}{2}(s_{m+1}c_{m+1}+s_mc_{m1})$$
(56)
The quantity $`s_m`$ has a sharp maximum for $`N1`$ at $`m_0=\frac{3\sqrt{5}}{4}N`$ and $`s_{m_0}`$ is

$$s_{m_0}=\frac{N}{2}(\sqrt{5}1)\sqrt{\sqrt{5}2}$$
(57)
The ground state energy in the thermodynamic limit is defined as a lowest eigenvalue of Eqs.(56), which is

$$E_0=\frac{N}{4}h_xs_{m_0}=\frac{N}{4}0.30028Nh_x$$
(58)
In the frame of variational function approach (51) the ground state magnetizations $`S_n^z`$ and $`S_n^x`$ are

$`S_n^z`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{m_0}{2N}}=0.309`$
$`S_n^x`$ $`=`$ $`{\displaystyle \frac{s_{m_0}}{N}}=0.300028`$ (59)
The function (51) has a momentum $`q=0`$. To calculate the spectrum in the first order in $`h_x`$ it is necessary to choose the variational function of the type (51) with momentum $`q`$. We omit here rather cumbersome calculations, which shows that the mass gap corresponds to $`q=\pi `$ and

$$m=h_x$$
(60)
The numerical diagonalization of finite cyclic systems shows very rapid exponential convergence to the thermodynamic limit. The extrapolated values for the ground state energy and the mass gap (at $`q=\pi `$) for $`h_x1`$ are

$`{\displaystyle \frac{E_0}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{4}}0.3017(1)h_x`$
$`m`$ $`=`$ $`0.4841(1)h_x`$ (61)
Comparison of these results with Eq.(58) shows that the variational ground state energy differs from the ‘exact’ one within $`0.4\%`$. As follows from Eq.(61) the mass gap in the multicritical point opens linearly with $`h_x`$, which is correctly described by the variational approach (60).

## V The critical line

In general, the critical line at $`h_x0`$ and $`h_z0`$ can not be found exactly. To obtain it we used the DMRG technique white . We have performed DMRG calculations using both the infinite-size and the finite-size DMRG algorithms. We calculated the ground state energy $`E_0(N)`$ and two lowest excitations $`m_1(N)=E_1(N)E_0(N)`$ and $`m_2(N)=E_2(N)E_0(N)`$.

In order to check the accuracy of the DMRG method we compared the obtained results with the exact ones for the ITF model. We used the infinite-size algorithm and open boundary conditions. The dependence of the results on the number of retained states $`s`$ in the DMRG computation and on a number $`N_{RG}`$ of DMRG steps ($`N=2N_{RG}+2`$) has been investigated. We have found that the calculation with $`s=25`$ gives satisfactory accuracy up to $`N=300`$. For example, relative errors in the ground state energy and in the mass gap at $`h_x=\frac{1}{2}`$ are $`\frac{\mathrm{\Delta }E_0}{E_0}=10^9(10^7)`$ and $`\frac{\mathrm{\Delta }m}{m}=10^7(10^5)`$ for $`N=100(300)`$. The accuracy becomes better when the value $`\left|h_x\frac{1}{2}\right|`$ is increased.

The critical field $`h_{x\mathrm{c}}(h_z)`$ at a fixed value $`h_z`$ ($`0<h_z<1`$) is determined by vanishing of the gaps $`m_1`$ and $`m_2`$. Below the critical field the mass gap $`m_1(N)0`$ exponentially with $`N`$. This behavior confirms the fact that the ground state is doubly degenerate in the thermodynamic limit at $`h_x<h_{x\mathrm{c}}(h_z)`$. The true mass gap in this region is defined by the value of $`m_2`$. The typical behavior of the gaps $`m_1`$and $`m_2`$ extrapolated to the thermodynamic limit is shown on Fig. 3. We note that at $`h_x>h_{x\mathrm{c}}(h_z)`$ the gaps $`m_1`$and $`m_2`$ coincide at $`N\mathrm{}`$.

The critical line obtained from the DMRG results extrapolated to the thermodynamic limit is shown on Fig. 1. In the plane $`(h_x,h_z)`$ it connects two limiting critical points studied in previous sections: the critical point ($`h_x=\frac{1}{2},h_z=0`$) of the ITF model and the multicritical point ($`h_x=0,h_z=1`$). Now we discuss the properties of the model (2) near the critical line.

The ITF model in the critical point is described by a conformal field theory with a central charge $`c=\frac{1}{2}`$. We expect that the ITF model is generic for the model (2) on the whole critical line (except the multicritical point $`h_x=0,h_z=1`$). To verify this suggestion we estimated the value of $`c`$ on the critical line. For this we used well known fact cardy  that for conformal invariant model with periodic boundary conditions the central charge appears at $`1/N`$ correction to the ground state energy

$$E_0=e_{\mathrm{}}N\frac{\pi c\upsilon }{6N}$$
(62)
where $`e_{\mathrm{}}`$ is the ground state energy per site in the thermodynamic limit and $`\upsilon `$ is the sound velocity.

At first, we calculated the sound velocity $`\upsilon `$ as (similar to the spectrum of the ITF model (18) the sound velocity on the whole critical line is determined at $`k=\pi `$)

$$\upsilon (N)=\frac{N}{2\pi }\left[E(\pi \frac{2\pi }{N})E(\pi )\right]$$
(63)
We carried out these calculations using the numerical diagonalization of the Hamiltonian (2) with the periodic boundary condition for $`N14`$. The size extrapolation is carried out by formula $`\upsilon (N)=\upsilon +aN^2`$. The example of the extrapolation procedure is shown on Fig. 4. The dependence of $`\upsilon `$ on $`h_z`$ along the critical line is shown on Fig. 5. After that, the central charge $`c`$ has been calculated using Eq.(62). The size extrapolation shows that the central charge is $`c=0.500(1)`$ for all calculated critical points.

The critical exponents $`\nu `$ of the ground state correlation function $`\widehat{O}_1\widehat{O}_RR^\nu (R1)`$ is defined by the scaling dimension $`X`$ of the operator $`\widehat{O}`$ and $`\nu =2X`$ book . The scaling dimension is related to the finite-size correction of the lowest excited eigenstate which can be reached from the ground state by the operator $`\widehat{O}`$. The finite-size correction to this excitation energy is

$$\mathrm{\Delta }E_i=E_iE_0=\frac{2\pi \upsilon }{N}X_i+o(\frac{1}{N})$$
(64)
For the ITF model at $`h_x=\frac{1}{2}`$ the sound velocity is $`\upsilon =\frac{1}{2}`$. The ground state has a momentum $`q=0`$ and the lowest excitation energy in the sector with momentum $`q=0`$ is

$$E(q=0)E_0=2\mathrm{sin}\frac{\pi }{2N}$$
(65)
and in the sector with $`q=\pi `$ is

$$E(q=\pi )E_0=\frac{1}{2}\mathrm{tan}\frac{\pi }{4N}$$
(66)
Thus, the scaling dimensions are $`X_0=1`$ and $`X_\pi =\frac{1}{8}`$. The corresponding associated operators are

$`\widehat{O}_{q=0}`$ $`=`$ $`{\displaystyle \underset{n}{}}S_n^x,`$
$`\widehat{O}_{q=\pi }`$ $`=`$ $`{\displaystyle \underset{n}{}}(1)^nS_n^z`$ (67)
which is in accord with the well known results McCoy  for the asymptotic of the correlation functions

$`S_0^xS_R^xS_0^x^2`$ $``$ $`R^2`$
$`S_0^zS_R^z`$ $``$ $`{\displaystyle \frac{(1)^R}{R^{1/4}}}`$ (68)
Using Eq.(64) we numerically checked that the scaling dimensions related to the finite-size corrections of the lowest excited eigenstates on the critical line remain as in ITF model $`X_0=1`$ and $`X_\pi =\frac{1}{8}`$. Therefore, we conclude that the model (2) on the critical line belongs to the universality class of the ITF model. This means that in accord with the prediction of the classical approach to the left of the critical line the staggered magnetizations $`(1)^nS_n^z`$ and $`(1)^nS_n^x`$ exist. But in contrast to the classical approach they vanish on the critical line with the critical exponent $`\frac{1}{8}`$. The mass gap is closed on the critical line and the critical exponent for the gap is equal to unity, i.e. the scaling behavior of the gap near the critical line is linear. For the particular case $`h_z=0.5`$, this is illustrated on Fig. 6. where the scaling plot for the scaled mass $`Nm_1`$ with $`N(h_xh_{x\mathrm{c}})`$ is shown.

The behavior of the critical line near the ITF point $`h_x=1/2,h_z=0`$ can be found from the following consideration. As it was established above in the vicinity of the critical line the gap is proportional to deviation from the line. This is valid for any direction of deviation except the direction at a tangent to the critical line. In the vicinity of the ITF critical point for fixed $`h_z1`$ according to Eq.(17) the gap is $`m=\left|h_xh_{x\mathrm{c}}(h_z)\right|`$. On the other hand, the gap is given by Eq.(36). Combining these two expressions for the gap on the line $`h_x=1/2`$, we obtain equation for the critical line in the vicinity of the point $`h_x=1/2,h_z=0`$ as

$$h_{x\mathrm{c}}(h_z)=\frac{1}{2}0.37h_z^2$$
(69)
This expression for the critical line is in a very good agreement with our numerical estimations up to $`h_z0.5`$.

In the vicinity of the multicritical point the form of the critical line can be found in a similar way. At $`h_z=1`$ and $`h_x1`$ the mass gap is proportional to $`h_x`$ (see Eq.(61)). On the other hand, for fixed $`h_x1`$ this gap is proportional to a deviation from the critical line $`m\left(h_zh_{z\mathrm{c}}\right)`$. Therefore, the critical line near the multicritical point behaves as $`h_{x\mathrm{c}}=A(1h_z)`$. The numerical coefficient $`A`$ has been found from the DMRG calculation of the critical line in the vicinity of the multicritical point. As a result, the form of the critical line at $`h_x1`$ is

$$h_{x\mathrm{c}}(h_z)=1.50\left(1h_z\right)+O\left(\left(1h_z\right)^2\right)$$
(70)
This form differs from that given in sen , where $`h_{x\mathrm{c}}(h_z)\left(1h_z\right)^{0.75}`$.

## VI Relation with the two-dimensional Ising model

It is well-known that the critical properties of the two-dimensional Ising model are equivalent to those of the one-dimensional ITF model kogut . Using the formalism of Ref.kogut  one can show that the quantum Hamiltonian (2) is related to the transfer matrix of the strong anisotropic (quasi-one-dimensional) Ising model in the uniform magnetic field. This model describes ferromagnetic Ising chains weakly antiferromagnetically coupled with each other. The Hamiltonian of the model is

$$H=J_1\sigma _{n,m}\sigma _{n+1,m}+J_2\sigma _{n,m}\sigma _{n,m+1}h\sigma _{n,m}$$
(71)
where $`\sigma _{n,m}=\pm 1`$, $`J_1J_2>0`$.

The relation between the model (2) and the model (71) are given by

$`2h_x\beta J_2`$ $`=`$ $`e^{2\beta J_1}1`$
$`h`$ $`=`$ $`2h_zJ_2`$ (72)
with $`\beta =1/kT`$.

The equivalence of these two models means that the model (71) undergoes the phase transition from the ordered (at $`T<T_c`$) to the disordered (at $`T>T_c`$) phase. The critical line $`T_c(h)`$ of the quasi-one dimensional model are related to the critical line $`h_{x\mathrm{c}}(h_z)`$ of the model (2). By use of mapping (72) we obtain that the critical temperature $`T_c(h)`$ is

$$kT_c(h)=\frac{2J_1}{\mathrm{log}\left(\frac{J_1}{J_2h_x(h/2J_2)}\right)}$$
(73)
In particular, at $`h0`$

$$kT_c=\frac{2J_1}{\mathrm{log}\left(\frac{2J_1}{J_2}\right)}$$
(74)
and $`T_c0`$ when $`h2J_2`$ as

$$kT_c=\frac{2J_1}{\mathrm{log}\left(\frac{J_1}{2J_2h}\right)}$$
(75)
The order parameter $`\frac{1}{2}\sigma _{n,m}+\sigma _{n,m+1}`$ vanishes at $`TT_c`$ with the critical exponent $`1/8`$. We note that the phase transition does not occur if $`J_2<0`$ as well as in the ferromagnetic version of the model (2).

The free energy of the model (71) is related to the ground state energy of the Hamiltonian (2) by

$$F=4J_2E_0(h_x,h_z)$$
(76)
According to Eqs.(49), (72) and (76) the zero-field susceptibility in the vicinity of the critical temperature $`T_c(0)`$ can be obtained using exactly solvable effective model (40). At $`\left|TT_c(0)\right|T_c(0)`$ the susceptibility is

$$\chi J_2=0.14120.24\frac{J_1}{kT_c(0)}\frac{TT_c(0)}{T_c(0)}\mathrm{log}\left(\frac{TT_c(0)}{T_c(0)}\right)$$
(77)
The last equation has a form coinciding with the results obtained for the two-dimensional antiferromagnetic Ising model by the series expansion method sykes .

## VII Summary

We have studied the antiferromagnetic Ising chain in the mixed transverse and longitudinal magnetic field. It was shown that the quantum phase transition existing in the ITF model remains in the presence of the uniform longitudinal field. Using the DMRG simulations we have found the critical line in the ($`h_x,h_z`$) plane where the mass gap is closed and the staggered magnetizations along the $`X`$ and $`Z`$ axes vanish. It is found numerically that the model on the critical line is described by the conformal field theory with the central charge $`c=1/2`$, i.e. it belongs to the universality class of the ITF model.

The scaling behavior in the vicinity of the ITF critical point is studied in detail. It is shown that the mass gap is proportional to $`h_z^2`$ and the contributions to it are given not only by the second order correction but also by all other divergent orders of the perturbation series in $`h_z`$. Besides, the analysis of the perturbation theory shows that the considered model at $`h_z1`$ can be mapped to the effective ITF model with renormalized parameters depending on $`h_z`$. In a framework of the effective ITF model the behavior of the susceptibility at $`h_x1/2`$ is determined.

The behavior of the model in the vicinity of the multicritical point is investigated. Using both the variational approach and the numerical diagonalization results we have found that the mass gap is proportional to $`h_x`$. Close to the multicritical point the form of the critical line is linear.

Of course, the considered model is the simplest case of the $`XXZ`$ model (1). It is interesting to extend the present analysis to this model as well as to study effects of interchain interactions in quasi-one-dimensional generalization of the model (2).

The mapping of the quantum model (2) to the strongly anisotropic statistical two-dimensional Ising model in the uniform magnetic field was considered. The behavior of the susceptibility of this model near the critical temperature is found. We expect also, that the main features of the considered model (2) are valid for the statistical two-dimensional antiferromagnetic Ising model in the uniform magnetic field $`h`$. That is, the applied magnetic field does not smear the phase transition existing in the two-dimensional Ising model, which is generic case for the whole transition line $`T_c(h)`$ in the plane ($`T,h`$). In particular, we believe that the analysis of the perturbation series in $`h1`$ in the same manner as was done in Sec.III gives the scaling behavior for the correlation length and the form of the critical line near the point $`T_c(0)`$ is $`T_c(0)T_c(h)h^2`$.

This work is supported under RFBR Grant No 03-03-32141 and ISTC No 2207.


format: text

text:
# Probing vortices in 4He nanodroplets

## Abstract

We present static and dynamical properties of linear vortices in 4He droplets obtained from Density Functional calculations. By comparing the adsorption properties of different atomic impurities embedded in pure droplets and in droplets where a quantized vortex has been created, we suggest that Ca atoms should be the dopant of choice to detect vortices by means of spectroscopic experiments.

The unique environment realized in liquid 4He clusters has opened up in recent years new opportunities for atomic/molecular spectroscopy to probe superfluid phenomena on the atomic scale toennies ; kwon . Helium droplets represent ideal nano-scale cryostats for a variety of fundamental experiments on liquid 4He, including the study of quantized vortices close . Vortices, while energetically unfavorable bauer , can potentially be stabilized by atomic or molecular impurities barranco1 .

During the free jet expansion experiments described in Refs. toennies ; close , it is plausible that quantized vortices may be created in some metastable state, long-lived enough to be detected. However, the question of whether 4He droplets can sustain vortices is still not resolved, and all the high resolution spectra of embedded molecules can apparently be explained without invoking their presence. Yet, it is expected that in the near future they could be created by some extension of the present experimental techniques. This calls for identifying signatures that might reveal vortical states in helium droplets. A possible experiment to detect their presence has been described by Close et al. close . They have suggested that alkali atoms, that normally reside in a ‘dimple’ on the surface of 4He clusters scoles ; anci1 ; ernst1 , may be drawn, when a vortex is present, inside the cluster along the vortex core. Spectroscopic experiments on the dopant atoms could thus provide evidence of their existence, since the line broadenings and shifts would be different in the two cases. We show in the following that alkali atoms are actually not suited for such an experiment, but rather alkaline earth (Ca) atoms may serve as probes to detect vortices.

Density Functional (DF) methods evans  have become increasingly popular in recent years as a useful computational tool to study the properties of classical and quantum inhomogeneous fluids, especially for large systems, for which they provide a good compromise between accuracy and computational cost. In particular, a quite accurate description of the properties of inhomogeneous liquid 4He at zero temperature has been obtained within a DF approach by using the energy functional proposed in Ref. dupont  and later improved in Ref. prica . This later DF, which has been successfully used over recent years to study a variety of 4He systems like clusters and films, is the one we use in the present work.

The minimization of the energy functional with respect to density variations, subject to the constraint of a given number of 4He atoms $`N`$, leads to the equilibrium particle density profile $`\rho (𝐫)`$, thus allowing to study the static properties of the 4He system. When dynamical properties are studied (as described in the following), we use the Time-Dependent DF (TDDF) method developed in Ref. giacomazzi , which allows to obtain both the 4He particle density $`\rho (𝐫,t)`$ and the velocity field $`𝐯(𝐫,t)`$. Briefly, in the static (dynamic) case, one has to solve a stationary (time-dependent) non-linear Schrödinger-like equation for an ‘order parameter’ $`\mathrm{\Psi }(𝐫)`$ $`(\mathrm{\Psi }(𝐫,t))`$, where the Hamiltonian operator is given by $`H=\frac{\mathrm{}}{2m}^2+U[\rho ,𝐯]`$. The effective potential $`U`$ is defined as the variational derivative of the energy functional, and its explicit expression is given in Ref. giacomazzi . From the knowledge of $`\mathrm{\Psi }\varphi e^{i\mathrm{\Theta }}`$ one can get the density $`\rho (𝐫,t)=\varphi ^2`$ and the fluid velocity field $`𝐯(𝐫,t)=\frac{\mathrm{}}{m}\mathrm{\Theta }`$. To model the interaction of liquid 4He with foreign impurities we use suitable He-impurity pair interaction potentials, which will be described later on.

We work in 3D cartesian coordinates, and adopt the following procedure to generate a vortex in the cluster in the most unbiased way. We consider a cluster in a rotating frame of reference with constant angular velocity $`\omega _z`$ around the z-axis don91 . The Hamiltonian density $`H`$ then acquires an additional term $`\omega _z\widehat{L}_z`$, $`\widehat{L}_z`$ being the angular momentum component along the z-axis. We minimize $`\mathrm{\Psi }`$ for this constrained Hamiltonian, imposing $`\mathrm{\Psi }`$ to be orthogonal, during the minimization, to $`\mathrm{\Psi }_0=\sqrt{\rho _{eq}(𝐫)}`$ describing the minimum energy state of the vortex-free cluster; we have applied the method to a $`N=300`$ droplet. To generate a vortex line, $`\omega _z`$ must be larger than a critical value -unknown in advance- $`\mathrm{\Omega }_c=\mathrm{\Delta }E/(N\mathrm{})`$ dal96 , where $`\mathrm{\Delta }E`$ is the energy cost to create a vortex (which in the present case is about 70 K, see Table 1, and hence $`\mathrm{\Omega }_c3\times 10^{10}`$ s-1), but not so large that one could generate a vortex array don91 . The particle density corresponding to this vortical configuration is shown in Fig. 1(a). We have calculated the circulation of the velocity field along a path enclosing the vortex core, and have exactly found the value $`Nh/m`$ appropriated for a quantized vortex with $`\widehat{L}_z=N\mathrm{}`$.

Note that since the vortex is quantized, the vortical state is an eigenstate of the angular momentum along the rotation axis, $`\widehat{L}_z`$. This means that our density profile is the same as that one would obtain by using the Feynman-Onsager ansatz, i.e. by adding to the energy functional an extra centrifugal term associated with an order parameter of the form $`\sqrt{\rho }e^{i\mathrm{\Phi }}`$ ($`\mathrm{\Phi }`$ being the azimuthal angle), and finding the density profile by solving an equation in the real quantity $`\sqrt{\rho (𝐫)}`$. This is the procedure used in the DF calculations of Ref. barranco1  to generate quantized vortex structures in helium drops -and also in Bose-Einstein condensates of trapped gases dal96 -. Instead, we have not assumed a priori a quantized value for the total angular momentum, but rather we generate a fully quantized vortex state starting from a pure cluster.

To use atomic impurities as probes of the presence of vortices in 4He drops, ideally one would like to have an atom that is barely stable on the surface of a pure drop, and becomes solvated in its interior in the presence of a vortex. The question of solvation vs. surface location for an impurity atom in liquid 4He can be addressed in an approximate way within the model of Ref. lerner  where, based on calculations of the energetics of impurities interacting with liquid 4He, a simple criterion is proposed to decide whether surface or solvated states are favored. An adimensional parameter is defined in terms of the impurity-He potential well depth $`ϵ`$ and the minimum position $`r_m`$, $`\lambda \rho ϵr_m/(2^{1/6}\sigma )`$, where $`\rho `$ and $`\sigma `$ are the bulk liquid density and surface tension of 4He, respectively. The criterion for solvation reads $`\lambda >1.9`$ for the existence of solvated states lerner . One thus needs an impurity with $`\lambda 2`$, and such that its most stable state is on the drop surface.

Alkali atoms are known to have their stable state on the surface of liquid 4He scoles ; anci1 , and they lie in the low $`\lambda `$ regime ($`\lambda 0.60.9`$) lerner . Accordingly, for alkalis a surface state should always be preferred, even in the presence of a vortex line. We have verified this point considering Na and Rb, as representative of light and heavy alkali, respectively. The alkali-He interaction is of the form proposed by Patil patil . We have compared the stable ‘dimple’ states of Na and Rb atoms on the surface of the $`N=300`$ cluster hosting a vortex line, with those of the same impurity trapped in the vortex core, exactly at the cluster center. We have found that the latter are energetically unfavored with respect to surface states, see Table 1. It is worth to note that in the case of Na, our results compare well with the Path Integral Monte Carlo calculations of Ref. nakayama , where a binding energy of about $`7`$ K is found for this cluster. We also note that, unlike the case of strongly bound impurities to 4He clusters barranco1 , which have their stable state inside the cluster, and for which there exists a critical cluster size below which the droplet+dopant+vortex complex is stable, the alkalis cannot stabilize the vortex, whatever the droplet size is barranco1 .

There are other dopants, however, for which there is clear evidence of a surface state on liquid 4He, i.e. alkaline earth atoms. Absorption spectra of alkaline earth atoms (Ca, Ba and Sr) attached to 4He clusters clearly support an outside location of Ca and Sr frank1 , and probably also of Ba frank2 . To describe the He-impurity interaction we employ an accurate ab-initio He-Ca pair potential meyer  used to study 4HeN+Ca droplets up to $`N=75`$ by Diffusion Monte Carlo techniques marius . For such a potential, $`\lambda 2.2`$, which apparently indicates a solvated stable state. However, for those cases where $`\lambda `$ is close to the solvation threshold, consideration of the shape of the potential energy surface, as well as the well depth and equilibrium internuclear distance, seems warranted scoles1 .

The stable state of a Ca atom in a 4He300 vortex-free cluster is shown in Fig. 1(b). Note that, in qualitative agreement with the experimental evidence, the ‘dimple’ appears to be much more pronounced than in the case of alkalis giacomazzi , reflecting the stronger He-atom interaction. However, in the presence of a vortex, the stable state is in the center of the cluster, as depicted in Fig. 1(c). The surface state in this case is unstable and, as the minimization proceeds, a Ca atom initially placed on the surface near the vortex core, is gradually drawn towards it and then sucked inside, eventually reaching the stable state in the center of the drop. The value of the angular momentum for the converged 4He configuration is again $`\widehat{L}_z=N\mathrm{}`$.

The response of the impurity atom to the different 4He environments shown in Fig. 1 might be determined with spectroscopic measurements, allowing to detect the presence of vortices: the observed linewidths and shifts of the excitation/emission spectra in the two cases shown in Fig. 1 should be very different, reflecting the ‘bubble’ environment in one case \[Fig. 1(c)\], and a more open environment in the other case \[Fig. 1(b)\]. Moreover, in the case of Fig. 1(b), bound-unbound transitions should be observed with a significant probability, thus implying a strong asymmetry in the observed spectra.

In the above picture one is assuming that the vortex is long-lived enough to allow a Ca atom, picked up randomly by the cluster, to diffuse close to the top of the vortex core and then to be drawn inside. We have no direct proof of the stability of the cluster+vortex complex on experimental time scales. However, we have indications that the cluster+vortex+dopant complex should be stable at least on the nanosecond time-scale. This conclusion comes from very long computing time simulations, using the TDDF method giacomazzi , to study the dynamics of the cluster+vortex+impurity complex. During these simulations, the impurity atoms were allowed to oscillate inside the ‘bubble’ in the cluster center \[see Fig.1(c)\], and the vortex line was always found to be stable, without showing any tendency to shrink, bend or migrate towards the surface of the cluster.

The solvated state for Ca in a vortex-free cluster is a stationary but unstable configuration against any displacement of the atom off the cluster center, only a few K in energy above the stable, ‘dimple’ state (see Table 1). This is a consequence of the borderline value of $`\lambda `$ for this impurity, and implies that in a real experiment, a fraction of Ca atoms might be trapped inside the clusters for fairly long times, even in the absence of vortices. For these atoms, the spectroscopic signals would be similar to those coming from Ca atoms trapped in the vortex core, making it difficult to discriminate between the two cases. A line shape calculation scoles  using as an input the 4He density profiles around the impurity might help to distinguish between solvated states of Ca with and without vortex. Since the extension of the method of Ref. scoles  -which is applied there to the simpler case of the monoelectronic alkali atoms- is rather involved for two-electron systems, we have not carried out such a calculation. We instead suggest additional measurements which may help to discriminate between the states shown in Fig. 1.

It appears from our calculations that the energy of the impurity-cluster system is rather insensitive to the location of Ca along the vortex core, once the atom is embedded in it. Consequently, vibrational modes of the impurity along the vortex line are expected to be soft. We have confirmed this by TDDF calculations, applying to the Ca atom a small initial momentum in a given direction -radially, towards the surface of the cluster for the ‘dimple’ state of Fig. 1(b), and along or perpendicular to the vortex line in the case shown in Fig. 1(c)-. We then let the impurity evolve in time, allowing for the 4He environment to dynamically follow the atom motion while the impurity oscillates around its equilibrium position. This is done in practice by numerically solving by means of a discrete Verlet algorithm, as usually done in Molecular Dynamics calculations, Newton’s equation of motion for the Ca atom under the force due to the surrounding 4He liquid

$$M\frac{d^2𝐑}{dt^2}=_𝐑\left\{\rho (𝐫,t)V_{HeCa}(𝐑(t)𝐫)𝑑𝐫\right\}$$

In this expression $`V_{HeCa}`$ is the pair potential describing the He-impurity interaction, and the density $`\rho (𝐫,t)`$ is updated at each time step according to the TDDF scheme for 4He giacomazzi . From the positions of the Ca atom as a function of time, relative to the center-of-mass of the Ca-He droplet system, different frequencies characterizing the impurity dynamics can be found from a Fourier analysis of the calculated time series.

We report in Fig. 2 the calculated vibrational spectra. The intensities are in arbitrary units, and normalized so that the higher peak in each spectrum has unit height note1 . It appears that the oscillation of Ca along the vortex core is indeed characterized by a single low frequency mode, as compared with the more fragmented spectrum for vibrations perpendicular to the vortex core. The presence of the soft mode is a signature of solvation of a Ca atom inside the vortex. Indeed, such a mode should be severely damped, or even absent, for a solvated Ca atom in a vortex-free cluster, since this configuration is unstable. We also show for comparison the vibration spectrum of a Ca atom in the ‘dimple’ state on the surface of a cluster without vortex. The peak just below 1 K is due to the ‘dipolar’ vibration of the impurity inside the semi-spherical(spherical) cavity in which it is trapped in the ‘dimple’(‘bubble’) state. Additional peaks appear in the ‘dimple’ spectra because of the coupling of the Ca motion with the surface modes of the 4He nanodroplet, which have similar frequencies (for instance, the lowest energy, $`l=2`$ quadrupolar mode of a pure 4He300 droplet occurs at $`0.6`$ K giacomazzi ). All these modes lie in the microwave frequency regime and there are experimental ideas to measure the corresponding vibrational frequencies ernst .

It is worth to see how these modes are coupled, which is particularly apparent when we displace the impurity perpendicular to the vortex core (dashed line in Fig. 2). To trigger this oscillation, we have given to the the impurity a kinetic energy of about 2 K, three times as much as in the other two cases. One may see that the spectrum displays one peak corresponding to the ‘dipolar’ mode discussed previously, and also softer modes of characteristics similar to those found in the other two cases. This coupling is possible because the time evolution is adiabatic.

Finally, we would like to emphasize that the dynamical behavior of He-impurity systems depends in a sensitive way on the details of the He-atom pair interaction. This calls for improving the available interaction potentials to strengthen the scenario described here, or to help in finding other atomic/molecular impurities which may serve as probes of the presence of vortices in 4He droplets.

We thank M.W. Cole, F. Dalfovo, W. Ernst, S. Hernández, G. Scoles, F. Stienkemeier, and F. Toigo for useful comments and discussions. We thank Prof. W. Meyer for the permission of using his unpublished Ca-He potential. F.A. acknowledges funding from MIUR-COFIN 2001 (Italy), and M.B. and M.P. from BFM2002-01868 and 2001SGR00064 (Spain).

Figure Captions

Figure 1: From left to right, for a 4He300 cluster, atomic equidensity lines corresponding to: (a) the vortical state; (b) the deep ‘dimple’ stable state of Ca, and (c) the solvated stable state of Ca in the presence of a quantized vortex. Each box side has 45 $`\mathrm{\AA }`$ lenght.

Figure 2: Calculated frequency spectra for the oscillations of a Ca atom in a 4He300 cluster. Solid line: solvated state \[see Fig.1(c)\], Ca moving along the vortex core; dashed line: solvated state, Ca moving perpendicular to the vortex core; dotted line: surface state \[see Fig.1(b)\], Ca moving in the radial direction.


text:
# Extreme blazars

## 1 Introduction

TeV astronomy is now leaving its infancy and becoming adult. The few BL Lacs (see Fig. 1) detected by the existing Cherenkov telescopes will become hundreds when the promised tenfold increased sensitivity (and lower energy thresholds) of the new generation of telescopes will be reached. It is a new electromagnetic window which has been opened to scrutiny: the fact that it can be done on ground means that it is one order of magnitude (at least) less expensive than satellites, and this in turn means fast progress. There are at least two science pillars involving TeV astronomy: the first is of course the study of TeV emitting objects, and the violent physics involved, while the second is the cosmological issue of estimating the intergalactic infrared, optical (and in future UV, when the detecting threshold will reach the 10 GeV range) backgrounds, through the study of the absorption imprinted in the intrinsic spectrum of the TeV emitting sources through the photon–photon collision process (e.g. Stecker & De Jager 1997). For these issues it is now important to find good candidates for TeV emission within a range of redshifts. As discussed below, we think that the best candidates are among low power BL Lacs, characterized by a spectral energy distribution (SED) with the first peak (in a $`\nu F_\nu `$ plot) at X–ray frequencies. We call them extreme blazars. We argue that the best TeV emitters are extreme BL Lacs with relatively strong radio emission, which measures the amount of seed photons required for the inverse Compton scattering process.

## 2 Extreme blazars

Besides the phenomenological approach leading to the blazar sequence discussed in Fossati et al. (1998) and in Donato et al. (2001), to understand why blazars have their SED controlled by their bolometric luminosity we need to model their spectra, in order to derive the physical parameters of the emitting region such as the size, the value of the magnetic field, the particle density and the beaming factor. Ghisellini et al. (1998) considered the sample of blazars detected by EGRET with some spectral information in the EGRET band. We applied a one–zone (leptonic) model, including (besides SSC) the inverse Compton scattering process with photons produced externally to the jet. We found a clear correlation between the energy of the electrons emitting at the peaks of the SED, $`\gamma _{\mathrm{peak}}`$, and the value of the magnetic plus radiation energy density measured in the comoving frame, $`U`$. The correlation was of the kind $`\gamma _{\mathrm{peak}}U^{1/2}`$, i.e. the radiative cooling rate at $`\gamma _{\mathrm{peak}}`$, $`\dot{\gamma }\gamma _{\mathrm{peak}}^2U`$, is constant and the same for all sources. Since the more extreme BL Lacs were under–represented in this work, Ghisellini, Celotti & Costamante (2002) later extended the range of parameters by including sources with more extreme values of $`\gamma _{\mathrm{peak}}`$, corresponding to low power BL Lacs. This was possible through the BeppoSAX observations of extreme blazars, covering the range 0.1–100 keV (i.e. all these BL Lacs were detected in the hard X–ray band).

We found a new branch of the correlation at high $`\gamma _{\mathrm{peak}}`$, with $`\gamma _{\mathrm{peak}}U^1`$. We interpreted it as the effect of a finite timescale of injection, $`t_{\mathrm{inj}}`$: in powerful blazars all particles cool in $`t_{\mathrm{inj}}`$ (fast cooling regime), while in low power BL Lacs the radiative cooling is less severe (slow cooling regime), and only the high energy particles can cool in this timescales. If we take a snapshot of the spectrum at the end of the injection, i.e. at the maximum of the flare, we then have two behaviors: in the fast cooling regime $`\gamma _{\mathrm{peak}}`$ is determined by the low energy cutoff of the injected distribution, since it is at this energy that a break in the final distribution occurs. In the slow cooling regime, instead, $`\gamma _{\mathrm{peak}}`$ corresponds to particles whose cooling time equals the injection time, and we obtain $`\gamma _{\mathrm{peak}}U^1`$.

These results suggest a two–step acceleration process: first a phase of “pre–heating” determining the low energy cut–off of the injected distribution, and then a rapid acceleration leading to a non–thermal energy distribution. This two–step process can operate in all blazars: their different SED simply reflects the different degree of radiative cooling. The typical energies produced by the pre–heating (in the range $`\gamma _{\mathrm{min}}`$10–103) correspond to the balance of the heating and cooling rates: this gives $`\gamma _{\mathrm{min}}U^{0.5}`$. We have $`\gamma _{\mathrm{peak}}=\gamma _{\mathrm{min}}`$ in powerful blazars, where electrons of all energies cool in a time $`t_{\mathrm{inj}}`$. Instead, when the cooling is less severe, the particle distribution (after $`t_{\mathrm{inj}}`$) will have a break at $`\gamma _{\mathrm{peak}}>\gamma _{\mathrm{min}}`$, determined by $`t_{\mathrm{cool}}(\gamma _{\mathrm{peak}})=t_{\mathrm{inj}}`$, leading to $`\gamma _{\mathrm{peak}}U^1`$.

Fig. 2 also shows that the parameters corresponding to different states of specific TeV sources obey the general trend (see the “tracks” connecting different points belonging to the same source). This is true even when the observed synchrotron peak frequency increases with the observed luminosity (as during the 1997 flare of Mkn 501; Pian et al. 1998), contrary to the general sequence discussed by Fossati et al. (1998).

### 2.1 Jet power of TeV sources

The fact that the emission region producing most of the observed luminosity is well localized allows us to estimate the total power carried by the jet. In fact, through modeling, we derive the size, the magnetic field, the particle density and the bulk Lorentz factor of the emission region, and these are the quantities needed to derive the power carried in the form of bulk kinetic energy and in Poynting flux. Fig. 3 shows our results, and we can compare TeV sources with the other blazars. As can be seen, TeV sources have the same kinetic powers of other BL Lacs, but significantly less Poynting flux. Contrary to more powerful blazars, we do not find strong constraints, on energetic grounds, to limit the number of electron–positron pairs of their jets. This is due to the fact that the mean electron random Lorentz factor is very high in TeV sources, making the relativistic mass of electrons to be roughly equal to the proton mass: it makes a little difference then to have an electron–proton or a pure pair plasma (compare $`L_\mathrm{e}`$ with $`L_\mathrm{p}`$ in Fig. 3).

## 3 Candidates for strong TeV emission

To emit TeV photons we need TeV electrons, which emit, by the synchrotron process, in the X–ray band (for typical magnetic fields of 0.1–1 Gauss and beaming factors $`\delta 10`$). Therefore sources with a strong synchrotron X–ray flux are obviously good TeV candidates. On the other hand, to have a large TeV luminosity through synchrotron self Compton (SSC), we need also synchrotron IR–optical seed photons, and this leads to a trade-off, since sources with a prominent synchrotron peak at hard X–ray energies have relatively little IR–optical emission. Fig. 4 shows how BL Lac objects (belonging to several samples) are located in the X–ray – radio flux plane. Filled diamonds corresponds to the already TeV detected BL Lacs: all of them are contained in the rectangle limited by the solid lines. We can see that all the TeV BL Lacs are HBL (high frequency peak BL Lacs, Giommi & Padovani 1994) characterized by a broad band radio to X–ray spectral index $`\alpha _{\mathrm{RX}}<0.75`$ (lower dashed line), but all of them are relatively bright radio sources. We believe that this occurs because the radio flux is a measure of the amount of IR–optical emission, providing the seed photons for the scattering. It can also be seen that there are several other sources in the vicinity of the already detected TeV ones. We considered them as the best TeV candidates, and we plot in Fig. 5 the expected flux at three different energy thresholds according to two models: a simple SSC model and the phenomenological model discussed in Fossati et al. (1998) with the modifications introduced in Costamante & Ghisellini (2002).

## 4 Extreme variability

Besides the famous example of very rapid variability in the TeV band shown by Mkn 421 in 1996 (Gaidos et al. 1996) there have been other examples of nearly factor 2 flux variations in a timescale of $``$5–20 minutes: in Mkn 501 in X–rays (Catanese & Sambruna 2000), BL Lac in X–rays (Ravasio et al. 2002) and Mkn 421 again in the TeV band (Aharonian et al. 2002). This very rapid variability is however not common, as indicated by the power spectra and structure functions (e.g. Fossati et al. 2000; Tanihata et al. 2003) calculated for the same sources, which suggests a more relaxed typical timescale of order of a few hours. The ultra–fast variations must be produced by very compact regions of size $`10^{14}`$ cm, several times smaller than the “normal” emission regions, and one interesting possibility is that they may correspond to “explosive” events occurring within larger and relatively steadier emission regions. Just for illustration, consider two colliding shells (as in the internal shock scenario), dissipating energy throughout their volumes, but also triggering reconnection of the magnetic field lines in small sub–volumes, dissipating large amounts of energy in the form of relativistic particles. Besides their own synchrotron flux, these particles will also scatter the “ambient” photons produced by the larger region, enhancing the inverse Compton contribution. If this idea is true, then we expect that these fast variations of the inverse Compton flux occur on top of flares (produced by the shell–shell collision event), and not during quiescent states (Ghisellini & Tavecchio, in prep.).

## 5 Internal shock scenario

The complex phenomenology of blazars and of blazar jets can be framed in a coherent picture, in which the central engine does not work continuously, but is injecting power in the jet intermittently, through shells of matter moving at slightly different velocities. These shells catch up at some distance from the black hole, and collide dissipating part of their kinetic energy. This “internal shock” scenario is one of the leading ideas for the generation of the prompt emission of gamma–ray bursts, and can work even better for blazar jets. As in gamma–ray bursts, the typical variability timescales in this scenario is linked with the initial distance $`\mathrm{\Delta }R`$ between two consecutive blobs, which should be of order of the Schwarzchild radius (they collide at a distance $`\mathrm{\Gamma }^2\mathrm{\Delta }R`$, but the observed time variability is shortened by a factor $`\mathrm{\Gamma }^2`$ by the Doppler effect). This scenario is easy to visualize, it explains the main characteristics of blazar jets and it allows numerical estimates (e.g. Ghisellini 1999; Spada et al. 2001; Tanihata et al. 2003).


text:
# 1 Introduction

## 1 Introduction

The distance to Pleiades became a controversial issue ever since Hipparcos (ESA 1997, Perryman et al. 1997) provided distance measurement to this cluster (van Leeuwen and Ruiz 1997, Pinsonneault, M. H., et al. 1998), as the new distance modulus was 5.3 mag rather than the traditional value of 5.6 mag. This implies that either something is wrong with the commonly used main sequence fitting method, or systematic errors in Hipparcos parallaxes are comparable to random errors, i.e. much larger than advertised. Either solution is uncomfortable for one or another group of astronomers. Percival et al. (2003) present the most recent attempt to modify the main sequence fitting to agree with Hipparcos parallaxes, and they provide references to past attempts. The most serious challenge to Hipparcos parallax to Pleiades was made by Narayanan and Gould (1999), who combined Hipparcos proper motions with known radial velocities, and using the conventional moving cluster method obtained Pleiades distance modulus 5.6 mag.

Fortunately, there is a simple classical method to resolve the dispute. Atlas, one of the brightest stars in Pleiades, with $`V=3.6`$ is known to be a binary. Hipparcos:

http://astro.estec.esa.nl/Hipparcos/HIPcatalogueSearch.html

gives astrometric orbital period as 290.6598 days, based on a periodic motion of the light centroid. The binary has been resolved by the Palomar Testbed Interferometer and the Mark III Optical Interferometer (Pan et al. 1999). The astrometric orbit has a semi-major axis of $`12.9\pm 0.1`$ mas and eccentricity of $`0.246\pm 0.006`$. The magnitude difference between the two components was measured to be 1.7 mag in V band. In their abstract Pan et al. (1999) write: ‘This leads to a distance for Atlas which is significantly larger than the value of $`116\pm 3`$ pc determined by the Hipparcos for the Pleiades cluster’. Unfortunately, there was no good radial velocity curve available at the time.

The only radial velocity data readily available are those compiled by Abt et al. (1965), and they are plotted with the Hipparcos period in Fig. 1. Presumably, these refer to the brighter of the two components. Fortunately, modern digital detectors and advanced methods of data analysis (Ruciński 2002) should provide very accurate radial velocity curves for both components as the star is bright and the observed radial velocity amplitude is large. Combining astrometric orbit with spectroscopic orbit is the most robust, purely geometrical method to measure distances with high accuracy and no serious systematic errors. This is a classical method. It was used, for example, to measure accurate distance to Hyades prior to Hipparcos (e.g. Torres et al. 1997). A recent very impressive application of this method provided the most accurate distance to the Sgr $`\mathrm{A}^{}`$: $`8.0\pm 0.4`$ kpc (Eisenhauer et al. 2003).

It is not clear why the determination of accurate astrometric orbit of Atlas (Pan et al. 1999) was not followed with a 1% accurate determination of the distance to Pleiades. No matter what the outcome, this would resolve one of the more dramatic controversies in modern astrophysics, with important implications either for the main sequence fitting method, or for future accurate space astrometry.

Acknowledgments. It is a great pleasure to acknowledge discussions with Dr. L. Eyer and Dr. M. Konacki

## REFERENCES

* Abt, H. A., Barnes, R. C., Biggs, E. S., and Osmer, P. S. 1965, Astrophys. J., 142, 1604.
* Eisenhauer, F. et a. 2003, astro-ph/0306220, , .
* ESA 1997, Hipparcos and Tycho Catalogs, ESA SP-1200 (Nordvijk: ESA), Vol. 3, .
* Nrayanan, V. K., and Gould, A. 1999, Astrophys. J., 523, 328.
* Pan, X. P., Shao, M., and Kulkarni, S. 1999, AAS Meeting 195, 132.03, .
* Percival, S. M., Salaris, M., Kilkenny, D. 2003, Astron. Astrophys., 400, 541.
* Perryman, M. A. C. et al. 1997, Astron. Astrophys., 323, L49.
* Pinsonneault, M. H., et al. 1998, Astrophys. J., 504, 170.
* Ruciński, S. 2002, Astron. J., 124, 1746.
* Torres, G., Stefanik, R. P., and Latham, D. W. 1997, Astrophys. J., 474, 256.
* Van Leeuwen, F., and Ruiz, C. S. H. 1997, in ‘Hipparcos’, Venice, ed. B. Battrick and M. A. C. Perryman (Noordwijk: ESA), 689, .


text:
# Robust polarization-based quantum key distribution over collective-noise channel

## Abstract

We present two polarization-based protocols for quantum key distribution. The protocols encode key bits in noiseless subspaces or subsystems, and so can function over a quantum channel subjected to an arbitrary degree of collective noise, as occurs, for instance, due to rotation of polarizations in an optical fiber. These protocols can be implemented using only entangled photon-pair sources, single-photon rotations, and single-photon detectors. Thus, our proposals offer practical and realistic alternatives to existing schemes for quantum key distribution over optical fibers without resorting to interferometry or two-way quantum communication, thereby circumventing, respectively, the need for high precision timing and the threat of Trojan horse attacks.

Quantum key distribution (QKD), such as the BB84 protocol proposed by Bennett and Brassard in 1984, allows two parties (Alice and Bob) to generate an arbitrarily long random secret key provided that they initially share a short secret key and that they have access to a quantum channel BB84 . As opposed to classical key distribution, the secrecy of the generated key does not rely on computational assumptions but simply on the laws of physics: as long as quantum mechanics holds, the information available to an eavesdropper (Eve) can be made arbitrarily small.

Photons are obvious candidates for mediators of quantum information since they are fast, cheap, and interact weakly with the environment. Both free air and optical fiber based QKD have been realized experimentally; see GRTZ2002  and HNDP2002  for reviews. Any experimental implementation of QKD naturally has to deal with the issue of noise in the quantum channel, which substantially complicates the security of QKD, as Eve may attempt to disguise her eavesdropping as noise from another source. Standard security proofs deal with channel noise, including photon loss, and show that Eve acquires essentially no information provided the noise rate is not too high. Higher noise rates mandate lower key generation rates, and once it becomes too large, secure key generation is impossible.

Building a viable quantum cryptographic system therefore depends on ensuring that the noise rate is low. The degree of freedom used to encode the information can be the polarization of the photon, its phase, or some combination of both. Purely phase-based schemes have been realized experimentally TRT1993  but require complex interferometric setups, high precision timing, and stable low temperatures. Interferometry becomes even more challenging with multi-photon states because of the difficulty of keeping phase coherence between the photons. A scheme which escapes some of these limitations using a clever encoding of key bits was proposed recently WASST2003 .

Polarization-based schemes also come with a disadvantage as optical fibers rotate polarizations of transmitted photons, and the degree of rotation fluctuates over time. If left untreated, this would result in an unacceptably high error rate. A number of proposals have been made to handle this source of errors; we present a new solution which is in some ways superior. Singlet states $`|\mathrm{\Psi }^{}=\frac{1}{\sqrt{2}}(|01|10)`$, where $`\{|0,|1\}`$ is any basis of the qubit Hilbert space, have the property that they are unchanged under equal rotations on both qubits; this is the defining property of a noiseless subspace. (If one’s qubits are the polarization degrees of freedom of single photons, as we assume here, then $`|0`$ and $`|1`$ can be taken to denote, for instance, the vertical and horizontal polarization states.) We present two protocols that take advantage of this property to encode key bits in four- or three-photon states. These states should be experimentally realizable, and form simple examples of a noiseless subspace or subsystem, respectively (also called a decoherence-free subspace and subsystem) ZR1997 ; KLV2000 ; KBLW2001 .

Free-space QKD is largely immune to the problem of polarization rotation: the coupling between the photons and the molecules in the atmosphere can be absorbed in a dielectric constant to very good approximation. The same unfortunately cannot be said about optical fiber. Rather, the dielectric constant acquires a spatial and temporal dependence, yielding an overall time dependent unitary transformation of the polarization state of a single photon, $`U(t)`$, as the net effect of the fiber. This varies on the time-scale of thermal and mechanical fluctuations of the fiber, the shortest of which we will refer to as $`\tau _{fluc}`$. If the time delay between the photons is small compared to $`\tau _{fluc}`$, the effect of this noise on the state of $`N`$ photons is well approximated by

$$\rho _N[U(t)]^N\rho _N[U(t)^{}]^N,$$
(1)
where $`t`$ now denotes the time of transmission. This is known as the unitary collective noise model ZR1997 .

There are several ways to deal with collective noise. The most obvious way is to continuously estimate the transformation $`U(t)`$ and systematically compensate for it. However, this requires an interruption of the transmission and, if the fluctuations become too rapid, the communication channel becomes useless. A second possibility — a phase-polarization hybrid which has been used successfully to realize QKD over $`67km`$ SGGRZ2002  — uses the Faraday orthoconjugation effect Martinelli1989  to autocompensate the effect of $`U(t)`$. Roughly speaking, if the transformation on the photon during its transmission from Bob to Alice is described by $`U(t)`$, it can be $`U(t)^1`$ when the photon is transmitted back from Alice to Bob, yielding no net transformation overall. The quantum information can be encoded by an extra phase transformation performed by Alice before returning the photon to Bob. Obviously, this technique only works if $`U(t)`$ is roughly constant throughout the transmission of the photon; this sets an upper limit of $`c\tau _{fluc}`$ to the distance over which QKD can be implemented with this scheme. 111However, with today’s technology, photon loss is a much more serious limitation to the distance over which QKD can be achieved.

Although such two-way quantum communication can eliminate collective noise, it allows for new attacks not possible against BB84. Since Alice receives and emits signals, it is possible for Eve to probe her laboratory — a technique known as the Trojan horse attack. There are many ways in which she can do this. She could add a weak signal to the channel at a slightly different frequency and recover some information about Alice’s phase transformation by subsequently filtering the output signal. Eve could also try to entangle an ancilla system with the signal before it enters Alice’s lab and perform a joint measurement on the two after Alice has retransmitted the signal. She could also intercept the signal and send a different signal to Alice, and thereafter measure the output to estimate the applied phase transformation, etc. Technical solutions for some of these attacks have been proposed. However, Eve has an enormous variety of attacks available, and to prove true information-theoretic security, one must assume that Eve has arbitrary technological power; for instance, she can outperform the best frequency filtering available to Alice and Bob. Because of its inherent use of two-way quantum communication, the protocol is formally quite different from the standard BB84 protocol, and proving its security may be quite difficult.222Note, however, that QKD must make use of two-way classical communication between Alice and Bob, and this presents no particular barrier to security proofs. Indeed, taking full advantage of two-way classical communication results in a protocol with substantially greater tolerance for noise GL2003 .

The schemes we propose here are purely polarization-based and cope with collective noise without resorting to two-way quantum communication. In the first protocol, the quantum information is encoded in a noiseless subspace, while in the second protocol, it is encoded on a noiseless subsystem.333Such encodings also obviate the need for Alice and Bob to share a reference frame (such as, for instance, a known relative alignment of their linear polarizers), as has been pointed out in the context of quantum communication in Ref. BRS2003 . A noiseless subspace is invariant under the action of the collective noise operation (here $`U^N`$). Any state within it is therefore unaffected (modulo a global unphysical phase) by the noise. When such a subspace does not exist, it may still be possible to find a set of density operators which are invariant under the effect of noise. These density operators instead form a noiseless subsystem on which pure states can be encoded.

The protocols we present use singlet states as building blocks. Information is encoded in the pairing of the photons; the various ways of organizing three or more photons into pairs provide us with different states with which to encode information. The photons must be distinguishable if different pairings are to correspond to different physical states: they need to be labelled in some way. Physically, this means that the photons must differ with respect to some degree of freedom. Here, the photons are assumed to differ in their time of arrival; they are spatially separated in the optical fiber. Furthermore, each bit is encoded on multi-photon states which must also be distinguished even in the presence of noise. Therefore, the multiplets of photons must be spatially separated by a distance greater than the separation between individual photons inside a multiplet. The fluctuation time $`\tau _{fluc}`$ needs only to be large with respect to the difference in the arrival times of the first and last photon of a multiplet.

It is crucial that information about the pairing resides only in the polarization state of the photon. For example, variations in the frequency of the photons can reveal pairing information: the frequencies of the two photons in each singlet must add up to the frequency of the pump, but the frequencies of photons in different singlets need not match. By measuring the energy of the photons, but not their polarization, Eve can learn about how they are paired without affecting the outcomes of Bob’s measurements. This energy signature of the pairing can be eliminated by filtering the photons before they leave Alice’s lab. Indeed, if the band-width of the filter is smaller than the band-width of the pump laser, almost no information about the pairing can be recovered by Eve. Similarly, we must wash out the phase relation between photon pairs. That is, the information must only be encoded on the relative order of the photons, not in their absolute time of arrival. This can be achieved by choosing the time delay between the photons in each multiplet at random from a certain range.

We will now present our two protocols: the first one encodes the quantum information on four-photon states while the second only requires three-photon states. We do not prove their security here, but rather point only to their similarities with a protocol proposed by Bennett in 1992 (B92) Bennett1992  which is known to be secure TKI2002 . We hope to provide a complete security proof in a later paper. The first protocol requires the definition of three normalized states of a photon quartet:

(2)
All these states correspond to pairs of singlet states: in $`|\psi _1`$, photons one and two form a singlet state and so do photons three and four. The two other states correspond to the two other possible ways of pairing four photons as is illustrated by the diagrams. The states are invariant under uniform rotations, so in any basis, these states can be decomposed into the superpositions noted above, where

$`|a`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0101+|1010)`$
$`|b`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0110+|1001)`$
$`|c`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(|0011+|1100).`$
It is straightforward to verify that $`|\psi _i|\psi _j|=1/2`$ for $`ij`$. It is therefore impossible to reliably distinguish any pair of these states, but it is possible to make a measurement that provides some information. Measuring the polarization of all four photons allows Bob to distinguish states $`|a`$, $`|b`$, and $`|c`$. Therefore, if Alice restricts her transmission to one of a pair of states, Bob can tell which of the two states she sent $`50\%`$ of the time. For example, suppose Alice transmits one of the pair $`\{\psi _1,\psi _2\}`$. When Bob measures either “0101” or “1010”, he can conclude that she sent $`|\psi _1`$. When he gets either “1100” or “0011”, he concludes she sent $`|\psi _2`$. Given any other outcome, Bob can not deduce with certainty which state she sent.

Note that Bob’s measurement constitutes a probabilistic error-free discrimination procedure Ivanovic  for any pair of the three states Alice might submit. That is, it is a measurement that with some probability identifies, without error, which of the pair was submitted, and with some probability yields an inconclusive outcome.

We now present the protocol.

Protocol 1

1. Alice chooses a random $`(4+\delta )n`$ bit string $`X`$ and a random $`(4+\delta )n`$ trit string $`B`$.
2. Alice encodes each bit $`\{0,1\}`$ of $`X`$ according to $`\{\psi _1,\psi _2\}`$ if the corresponding trit of $`B`$ is $`0`$; $`\{\psi _2,\psi _3\}`$ if $`B`$ is 1; or $`\{\psi _3,\psi _1\}`$ if $`B`$ is 2.
3. Alice sends the $`(4+\delta )n`$ quartets of photons to Bob.
4. Bob receives the photons, and announces this fact. For each of the $`(4+\delta )n`$ photon quartets, he randomly chooses between the rectilinear or diagonal polarization basis. He then measures each of the four photons of each quartet according to this choice of basis.
5. Alice announces B. Given this information, and using the procedure described above, Bob can determine, for each quartet, whether or not his measurement was conclusive, and if conclusive, the value of the encoded bit.
6. Alice and Bob discard all bits where Bob’s measurement was inconclusive. With high probability, there are at least $`2n`$ bits left which they keep. Otherwise, they abort the protocol.
7. Alice selects a random subset of $`n`$ bits and tells Bob which bits were selected.
8. Alice and Bob announce and compare the value of the $`n`$ selected bits to estimate Eve’s interference; if more than an acceptable number of errors are found, they abort the protocol.
9. Alice and Bob perform information reconciliation and privacy amplification on the remaining $`n`$ bits.

In step 4, the choice of basis does not affect the measurement outcome of Bob: this is in fact the main property of the encoding. Nevertheless, it is crucial that Eve does not know in which basis the measurement is performed. If she knew, she could measure in the same basis as Bob, and would know everything Bob knew. Since she does not know, she will frequently measure in a different basis than Bob and therefore introduce errors that will reveal her presence.

As the protocol is written, in step 6, Alice and Bob discard any bits for which Bob’s measurement is inconclusive. Nevertheless, an inconclusive result could still be useful. Indeed, any measurement result whose weight differs from $`2`$, e.g. “1011”, indicates that Eve has tampered with the communication. This provides Alice and Bob with some extra data to estimate Eve’s interference: only allowed code-words should be observed by Bob.

If steps 4 and 5 are inverted — which could only provide Eve with more information — we get a protocol quite similar to B92. Alice encodes the bit in two preselected nonorthogonal states which she sends down the quantum channel. Bob then performs a von Neumann measurement chosen at random from a certain set, which can be cast in terms of a positive operator valued measurement (POVM): this is also required in B92. Nevertheless, there are certain differences in the nature of these POVMs which must be studied carefully to arrive at a complete security proof. We are currently working on these issues. The B92 protocol is not secure if the transmission rate is below 1/2. Eve can replace the noisy channel by a perfect channel and measure Alice’s output in such a way that she achieves a conclusive discrimination with probability , in which case, she knows the state and can send it to Bob. In the case of an inconclusive result, she does not send anything. From Alice and Bob’s point of view, this would be indistinguishable from the natural noise. By delaying the announcement of $`B`$, our proposal escapes this limitation.

The simplicity of the measurement — single photon polarization — is a clear advantage of this protocol. (The possibility of discriminating between states that are invariant under collective noise, without having recourse to collective measurements, has also been noted in Ref. Cabello2003 ) Furthermore, the required states have already been produced by several groups PDGWZ2001 ; ZYCZZP2003 . Two singlet states can be produced in a short time interval via parametric down-conversion by sending a femto-second pump laser pulse back and forth across a crystal (using a mirror). Since the photons are emitted in different directions, the EPR pairs can be clearly distinguished. Zhao et al. ZYCZZP2003  have produced these pairs of singlets at a rate of 20 kHz. Optical delays and switches can be used to create any of the three states of Eq. 2. If Bob’s single-photon detectors have a jitter of roughly 300 picoseconds, the delay lines must have a length of order 10 centimeters.

At first glance, it looks like there are two copies of the information in this encoding. For instance, if Alice announces that the bit is encoded as $`\{\psi _1,\psi _2\}`$, the value of the first two measurement outcomes is enough to sometimes deduce the value of the encoded bit: it is necessarily 1 if the outcomes are the same. The same holds for the measurement outcomes on the last two photons. Nevertheless, this redundancy is intrinsic to our quantum encoding scheme and doesn’t provide Eve with any extra information. The second protocol we present exploits this redundancy to reduce the size of the encoding.

Protocol 2 is a slight modification of Protocol 1. In step 3, instead of sending the entire state to Bob, Alice randomly discards one photon from each quartet and sends the remaining three. In step 5, Alice should also announce which photon she has discarded. Therefore, the three pure states of Eq.2 are replaced by the three mixed states

(3)
where the “$``$” denotes the maximally mixed state. These states are obviously invariant under collective noise. Furthermore, any pair can be distinguished with a finite probability, just as with the states in Protocol 1. This follows from the fact that they constitute non-orthogonal mixed states with non-identical supports, and the fact that one can achieve probabilistic error-free discrimination of such mixed states RST2003 . For example, suppose that Alice has announced that the bit is encoded as $`\{\rho _1,\rho _2\}`$. Clearly, any measurement outcome of Bob’s where photon 1 and photon 2 come out parallel rules out the state $`\rho _1`$. Here, the two outcomes “000” and “111” never occur in the absence of eavesdropping; they indicate that Eve has tampered with the communication.

These states do not form noiseless subspaces, because any particular pure state in the decomposition of the density matrices $`\rho _1`$, $`\rho _2`$, or $`\rho _3`$ does not remain invariant under collective noise. Instead it is transformed into another state in the decomposition of the same density matrix. For instance, $`|\mathrm{\Psi }^{}|0`$, in the decomposition of $`\rho _1`$, becomes under collective bit-flip $`|\mathrm{\Psi }^{}|1`$. The individual states are not noiseless, but the space spanned by them is; therefore they form a noiseless subsystem, and the density matrices are invariant under collective noise.

This second protocol shares many similarities with Protocol 1. It is our hope that essentially the same proof will be able to show both protocols are secure. In practice, Alice does not have to create two photon pairs and discard one photon; she could simply create one pair and one additional photon in the maximally mixed state. This should greatly increase the optimal transmission rate since pair creation schemes have relatively low efficiency. Furthermore, Protocol 2 — based on trios of photons instead of quartets — should suffer less from photon loss and hence be realizable over greater distances.

On a speculative note, perhaps the two protocols could be hybridized into a more robust protocol. Alice could always encode her information on photon quartets. Bob could then divide the photon multiplets into two sets depending on how many photons from the quartet actually made it to the destination: set 1 when all four photons made it and set 2 when only 3 photons were detected. The outputs from set 1 could then be used to complete Protocol 1 while those of the second set would be used as in Protocol 2. However, this suggestion provides Eve with a wide variety of attacks unavailable in our two protocols, and its security must therefore be studied on a different basis.

We thank Gilles Brassard for helpful comments on our protocols. J.-C.B., R.L., D.P., and R.W.S. receive financial support from Canada’s NSERC. This work was partially funded by MITACS and ARDA.