text:
# Is there an imprint of Planck scale physics on inflationary cosmology?

## I Introduction

The present upper limits on the temperature fluctuations of the microwave background radiation, corresponding to the density perturbations at the time of recombination, are less than $`1`$ part in $`10^5`$ on scales larger than the Hubble radius. On these scales, the observations of the CMBR gives the power spectral index $`n`$ to be $`1.1\pm 0.1`$. These observations provide an upper bound on the amplitude of density fluctuations at that epoch. Inflationary paradigm provides a physical mechanism for the generation of seed density perturbations (see, for example, Ref. linde ). Many models of inflation require a period much longer than 70 e-foldings to solve the horizon and flatness problems of standard cosmology. In order to obtain the observed density inhomogeneities on the galactic scales, these models of inflation require the modes of vacuum fluctuations at the initial epoch of inflation to be smaller than Planck length $`L_P(G\mathrm{}/c^3)^{1/2}`$ paddyold .

It is well known that the space-time structure at Planck scales is affected by the quantum gravitational effects and it is generally believed that $`L_P`$ acts as a physical cutoff for space-time intervals. The existence of a fundamental length implies that processes involving energies higher than Planck energies will be suppressed, and the ultraviolet behavior of the theory will be improved. All existing models of quantum gravity provide a mechanism for good ultra-violet behavior, essentially through the existence of the fundamental length scale. Even though the models of quantum gravity have some success, none of these have given a complete theory that works at the Planck scales.

If the ultimate theory of quantum gravity has a fundamental length scale, then the low-energy effective quantum field theories should have an imprint of the Planck scale. However, the standard formulation of quantum field theory does not take into account the existence of any fundamental length in the space-time. Attempts have been made to incorporate quantum gravitational effects with the help of two different approaches: (i) by introducing the fundamental length scale in a Lorentz-invariant manner in the Feynman Green function and (ii) by modifying the dispersion relation of the linear field equation.

In this work, our focus is on approach (ii). (For method (i) and its implications in conventional quantum field theory, see Ref. paddy1 ; paddy2 .) In the context of inflationary cosmology, trans-Planckian effects have gained interest after the detailed work of Martin and Brandenberger branden ; branden1  who introduced modified dispersion relations in an adhoc manner. The authors dealt with high-frequency dispersion relations, used earlier by Unruh unruh  and Corley/Jacobson corley  for the analysis of deviations from the thermality of the Hawking spectrum, to the inflaton field and obtained the spectrum of the vacuum fluctuations after it crosses the Hubble radius (during the inflation). They found that while some hypothesized dispersion relation(s) to quantum field theory yield no change in the power spectrum of the fluctuations, some indicate the cosmological sensitivity to Planck scale physics. However, there have been arguments that the spectrum of vacuum fluctuations should not be sensitive to the Planck scale physics staro ; paren .

All previous analyses of trans-Planckian effects on inflationary perturbations branden -other\_pap  have concentrated on specific form of the dispersion relations for a particular choice of initial state of the field. \[Two sets of initial conditions which are used in the literature are Minimum energy and Instantaneous Minkowski vacuum.\] In our analysis, we do not assume any specific form of the dispersion relation and also allow initial state of the field to be arbitrary. The only condition on the dispersion relation being that the non-linear dispersion relations become linear for $`\lambda >>L_P`$. From the earlier analyses, we would expect that the spectrum of fluctuations will strongly depend on the nature of the dispersion relation and the initial state of the field.

Using this strategy, we address the following questions: Is there an imprint of Planck scale physics on inflationary cosmology? and Can the deviations from the flat scale-invariant spectrum due to non-linear dispersion relations be differentiated from that of deviations from the flat spectrum in the standard inflationary scenario? We try to provide an answer by considering a quantum field satisfying the linear wave equation with higher spatial derivative terms propagating in the de Sitter space-time. We perform a general analysis for a class of dispersion relations whose field modes can be expanded as WKB modes over a range of initial states (conditions) for the quantum field. \[In this paper we are primarily interested in frequency evolution of modes and we ignore the actual creation of the modes. There are claims in the literature that the energy density of fluctuations of the inflaton field is necessarily high and its back-reaction on the metric cannot be ignored tanaka ; staro  (see also, Ref. tpsingh ).\]

We show that the signatures of the deviations from the flat scale-invariant spectrum due to non-linear dispersion relations cannot be differentiated from the standard inflationary scenario with an arbitrary initial state. We supplement this claim by studying four dispersion relations:

$`\mathrm{\Omega }^2`$ $`=`$ $`{\displaystyle \frac{E_P^2}{2\pi }}\mathrm{log}\left(1+{\displaystyle \frac{2\pi k^2}{E_P^2}}\right)`$
$`(\mathrm{Padmanabhan}^{}\mathrm{s}\mathrm{dispersion}\mathrm{relation})`$
$`\mathrm{\Omega }^2`$ $`=`$ $`\kappa ^2\mathrm{log}^2(1+k/\kappa ),`$
$`(\kappa \mathrm{Poincare}\mathrm{dispersion}\mathrm{relation})`$
$`\mathrm{\Omega }^2`$ $`=`$ $`k_{\mathrm{Pl}}^2\mathrm{tanh}[(k/k_{\mathrm{Pl}})^n]^{(2/n)},`$
$`(\mathrm{Unruh}^{}\mathrm{s}\mathrm{dispersion}\mathrm{relation})`$
$`\mathrm{\Omega }^2`$ $`=`$ $`k_{\mathrm{Pl}}^2\mathrm{tan}^2[k/k_{\mathrm{Pl}}],`$
$`(\mathrm{Ahluwalia}^{}\mathrm{s}\mathrm{dispersion}\mathrm{relation})`$
where $`E_P`$ is the Planck energy and $`\kappa ,k_{\mathrm{Pl}}`$ is the inverse Planck length (or of the order of Planck scale). The dispersion relation (I) reproduces the exact density of states required for the black-hole by treating black-hole as a one-particle state of a non-local field theory paddy98 , relation (I) is motivated by $`\kappa `$-Poincare algebra gilkman  while the last relation is the gravitationally modified expression for the wave-particle duality aluwalia .

The organization of the paper is as follows: In Sec. (II), we calculate the power spectrum of the vacuum fluctuations for an arbitrary initial condition in (i) standard inflationary scenario (ii) dispersive field theory models, and interpret the results. Finally, in Sec. (III) we summarize the results.

## II Power spectrum of quantum fluctuations

In this section, we calculate the power spectrum of fluctuations for an arbitrary initial state of a quantum field $`\mathrm{\Phi }`$, in standard inflationary scenario and in dispersive field theory models, propagating in the spatially flat FRW background (specifically de Sitter). We show that using the observations of CMBR the signatures of the deviations from the flat scale-invariant spectrum due to non-linear dispersion relations cannot be differentiated from the standard inflationary scenario with an arbitrary initial state.

### II.1 Standard inflation

Let us consider scalar fields propagating in the flat FRW metric of the form

$$ds^2=dt^2a^2(t)d𝒙^2=a^2(\eta )\left[d\eta ^2d𝒙^2\right],$$
(5)
where $`𝒙`$ is the 3-space, $`t`$ is the comoving time, $`\eta `$ is the conformal time and $`t=a(\eta )𝑑\eta `$. For flat de Sitter $`a(\eta )=1/(H\eta )`$ \[$`a(t)=\mathrm{exp}(Ht)`$\] and $`H`$ is the Hubble constant during inflation. \[Typical values of the inflationary expansion are $`E10^{15}GeV`$, $`H5\times 10^{24}cm^1`$, start of the inflation $`\eta _i10^{36}sec`$ ($`t_i10^{36}sec`$) and the end of inflation $`t_f70H^1`$.\] The free massless, minimally coupled scalar field satisfies the Klein-Gordon equation,

$$\mathrm{}\mathrm{\Phi }\frac{d^2\mathrm{\Phi }}{d\eta ^2}+2\frac{1}{a^3(\eta )}\frac{da(\eta )}{d\eta }^2\mathrm{\Phi }=0.$$
(6)
The symmetry of the Robertson-Walker metric allows for separating variables in Eq. (6) and the scalar field can be decomposed as

$$\mathrm{\Phi }(𝒙,\eta )=\frac{1}{(2\pi )^{3/2}}_k\left[a_ku_k(𝒙,\eta )+a_k^{}u_k^{}(𝒙,\eta )\right]𝑑𝒌,$$
(7)
where

$$u_k(𝒙,\eta )=a(\eta )^1\chi _k(\eta )\mathrm{exp}(i𝒌.𝒙).$$
(8)
\[$`𝒌`$ is the comoving wave vector and $`k^2=|𝒌|^2`$.\] Note that the creation and annihilation operators $`a_k^{},a_k`$ obey the usual commutation relations. The time-dependent part of the mode function satisfies the oscillatory equation

$$\frac{d^2\chi _k}{d\eta ^2}+\varpi _k^2(\eta )\chi _k(\eta )=0,$$
(9)
where $`\varpi _k^2(\eta )=k^2[a^{\prime \prime }(\eta )/a(\eta )]`$. \[The solution to the above equation reduces to the usual harmonic oscillator solution, $`\chi _k\mathrm{exp}(i\omega \eta )`$, with a linear dispersion relation, $`\omega ^2=k_{\mathrm{ph}}^2=k^2/a^2(\eta )`$, on scales much smaller than the Hubble length. $`k_{\mathrm{ph}}`$ is the physical momentum of the perturbation.\] A complete set of mode solutions to Eq. (9) is obtained by imposing initial conditions $`\chi _k(\eta _i)`$, $`\chi _k^{}(\eta _i)`$ on a Cauchy surface $`\eta =\eta _i`$ corresponding to a homogeneous initial state.

The solution to the above equation in the short-wavelength limit ($`\lambda _{\mathrm{ph}}<<H^1`$) corresponds to the Minkowski modes, i. e.,

$$\chi ^I=\frac{a_k}{\sqrt{2k}}\mathrm{exp}(ik\eta )+\frac{b_k}{\sqrt{2k}}\mathrm{exp}(ik\eta )$$
(10)
where, $`a_k`$ and $`b_k`$ are determined by the choice of the initial conditions. Orthonormality of the modes gives,

$$|a_k|^2|b_k|^2=1$$
(11)
Thus, $`a_k`$ and $`b_k`$ can be any arbitrary functions of $`k`$ satisfying the above condition. In the long wave-wavelength limit ($`\lambda _{\mathrm{ph}}>>H^1`$), the modes are frozen and the solution is given by

$$\chi ^{II}=C_ka(\eta )=\frac{C_k}{H\eta }.$$
(12)
Demanding that $`\chi ^I`$ and $`\chi ^{II}`$ are continuous at $`\eta _H`$, $`C_k`$ can be determined in-terms of $`a_k`$ and $`b_k`$. \[The junction condition at $`\eta _H`$ between the two modes correspond to the relation $`\{2\pi a(\eta )\}/k=H^1`$, i. e., $`\eta _H=2\pi /k`$.\] Thus, for the $`k`$ mode leaving the Hubble radius at $`\eta _H`$, we get

$$k^3|C_k|^2=\frac{H^2}{2}+H^2\left[|b_k|^2+|b_k|(|b_k|^2+1)^{1/2}\right],$$
(13)
where we have used condition the (11) to write $`|a_k|^2`$ in-terms of $`|b_k|^2`$. The above result shows that the spectrum of fluctuations has a scale-invariant term plus a non scale-invariant term which strongly depends on the choice of the initial state ($`b_k`$). If we assume that on cosmologically relevant scales the inflaton field is in the vacuum state, corresponding to no inflaton particles corresponding to the initial condition,

$$\chi (\eta _i)=\frac{1}{\sqrt{2k}};\chi ^{}(\eta _i)=\pm i\sqrt{\frac{k}{2}},$$
(14)
we get $`a_k=1,b_k=0`$; which results in a flat scale-invariant spectrum.

### II.2 Dispersive field theory model

As mentioned in the introduction, in many models of inflation, the physical length of the perturbations is smaller than the Planck length at the initial epoch implying we need to consider quantum gravitational corrections. One way of introducing quantum gravitational corrections is to modify the linear dispersion relation in an ad-hoc manner by breaking the Lorentz invariance i. e., by rewriting $`\omega ^2=\mathrm{\Omega }^2[k_{\mathrm{ph}}(\eta )]`$. Replacing $`k^2`$ in Eq. (9) as

$$k_{\mathrm{mod}}^2=a^2(\eta )\mathrm{\Omega }^2[k_{\mathrm{ph}}(\eta )]=a^2(\eta )\mathrm{\Omega }^2\left(\frac{k}{a(\eta )}\right),$$
(15)
we get,

$$\frac{d^2\chi _k}{d\eta ^2}+\varpi _{}^{\mathrm{mod}}{}_{k}{}^{2}(\eta )\chi _k=0,$$
(16)
where,

$$\varpi _{}^{\mathrm{mod}}{}_{k}{}^{2}(\eta )=a^2(\eta )\mathrm{\Omega }^2\left(\frac{k}{a(\eta )}\right)\frac{a^{\prime \prime }(\eta )}{a(\eta )}.$$
(17)
The modified dispersion relations must satisfy the property $`k_{\mathrm{mod}}k`$ for $`k<<k_{\mathrm{Pl}}`$. We have tabulated the three regimes of interest and the dominant contribution of $`\varpi _k^{\mathrm{mod}}`$ in these regimes in Table 1. Using the notation in Table I, the solution to Eq. (16) in regime I, in the WKB limit is

$`\chi _k^I(\eta )`$ $`=`$ $`{\displaystyle \frac{A_k}{\sqrt{2W_k^I(\eta )}}}\mathrm{exp}\left(i{\displaystyle W_k^I(\eta )𝑑\eta }\right)`$ (18)
$`+`$ $`{\displaystyle \frac{B_k}{\sqrt{2W_k^I(\eta )}}}\mathrm{exp}\left(i{\displaystyle W_k^I(\eta )𝑑\eta }\right),`$
where, $`A_k`$ and $`B_k`$ are constants and are determined by the choice of the initial conditions. Orthonormality of the modes gives,

$$|A_k|^2|B_k|^2=1.$$
(19)
Here again, we assume $`A_k`$ and $`B_k`$ to be arbitrary functions of $`k`$ satisfying the above condition. In regime II, solution to Eq. (16) is (Minkowski) plane wave solution, i. e.,

$$\chi _k^{II}(\eta )=\alpha _k\mathrm{exp}(ik\eta )+\beta _k\mathrm{exp}(ik\eta ).$$
(20)
The junction conditions of the wave modes $`\chi _k^I/\chi _k^{II}`$ and its derivatives at $`\eta =\eta _{\mathrm{Pl}}`$ gives

$`\alpha _k`$ $`=`$ $`A_k\vartheta _k\mathrm{exp}[iI(\eta _{Pl})+ik\eta _{Pl}],`$ (21)
$`\beta _k`$ $`=`$ $`B_k\vartheta _k\mathrm{exp}[iI(\eta _{Pl})ik\eta _{Pl}],`$ (22)
where,

$`\vartheta _k`$ $`=`$ $`\sqrt{{\displaystyle \frac{W_k^I(\eta _{Pl})}{8k^2}}}+{\displaystyle \frac{1}{\sqrt{8W_k^I(\eta _{Pl})}}},`$ (23)
$`I(\eta )`$ $`=`$ $`{\displaystyle _{\eta _i}^\eta }W_k^I(\eta ^{})𝑑\eta ^{}={\displaystyle _{\eta _i}^\eta }a(\eta ^{})\mathrm{\Omega }[k_{ph}(\eta ^{})]𝑑\eta ^{}.`$ (24)
In regime III, the solution is

$$\chi _k^{III}(\eta )=C_k^{\mathrm{mod}}a(\eta ),$$
(25)
where $`C_k^{\mathrm{mod}}`$ is a constant whose modulus square gives the power spectrum of the density perturbations. Using the junction conditions at $`\eta =\eta _H`$, we get

$$C_k^{\mathrm{mod}}=\frac{1}{a(\eta _H)}\left[\alpha _k\mathrm{exp}(ik\eta _H)+\beta _k\mathrm{exp}(ik\eta _H)\right].$$
(26)
Rewriting $`\alpha _k/\beta _k`$ in-terms of $`A_k`$ and $`B_k`$, we get

$`|C_k^{\mathrm{mod}}|^2=a^2\vartheta _k^2[|A_k|^2+|B_k|^2`$ (27)
$`+`$ $`(A_kB_k^{}\mathrm{exp}[2iI(\eta _{Pl})+2ik(\eta _{Pl}\eta _H)])]`$
Using the condition (19), we obtain

$`|C_k^{\mathrm{mod}}|^2={\displaystyle \frac{H^2k^2}{4k_{\mathrm{mod}}(\eta _H)}}(1+{\displaystyle \frac{k_{\mathrm{mod}}(\eta _H)}{k}})^2[{\displaystyle \frac{1}{2}}+|B_k|^2`$
$`+\sqrt{|B_k|^4+|B_k|^2}\mathrm{cos}[2I(\eta _{Pl})+2k(\eta _{Pl}\eta _H)]].`$ (28)
Using the property that $`k_{\mathrm{mod}}k`$ for $`k<<k_{\mathrm{pl}}`$, we get

$`k^3|C_k^{\mathrm{mod}}|^2`$ $`=`$ $`{\displaystyle \frac{H^2}{2}}+H^2|B_k|^2+H^2|B_k|\sqrt{|B_k|^2+1}`$ (29)
$`\times `$ $`\mathrm{cos}\left[2I(\eta _{Pl})+2k(\eta _{Pl}\eta _H)\right].`$
The above result shows that the spectrum of fluctuations has a scale-invariant term plus a non scale-invariant term which strongly depends on the initial state ($`B_k`$) and the form of the dispersion relation ($`\mathrm{\Omega }`$). If we assume the initial state of the field (Minimum energy state) to be

$`\chi _k(\eta _i)_{\mathrm{mod}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{2a(\eta _i)}}}\mathrm{\Omega }^{1/2}\left({\displaystyle \frac{k}{a(\eta _i)}}\right),`$ (30)
$`{\displaystyle \frac{1}{i}}\chi _k^{}(\eta _i)_{\mathrm{mod}}`$ $`=`$ $`\pm \sqrt{{\displaystyle \frac{a(\eta _i)}{2}}}\mathrm{\Omega }^{1/2}\left({\displaystyle \frac{k}{a(\eta _i)}}\right),`$ (31)
we get $`A_k=\mathrm{exp}[I(\eta _i)],B_k=0`$; which results in a flat scale-invariant spectrum.

### II.3 Analysis and Interpretation of results

In subsections (II.1) and (II.2), we obtained the power spectrum of fluctuations for an arbitrary initial state in standard inflationary scenario and dispersive field theory models. In this subsection, we compare and interpret the results obtained. We first show that the third term in the right hand side of (29) cancels out. We then show that the signatures of the deviations from the flat scale-invariant spectrum due to non-linear dispersion relations cannot be differentiated from the standard inflationary scenario with an arbitrary initial state.

Using the property that the non-linear dispersion relations becomes linear for $`k_{\mathrm{Pl}}>>k_{\mathrm{ph}}`$, we can write them in the general form

$$\mathrm{\Omega }=k_{\mathrm{Pl}}F\left(\frac{k_{\mathrm{ph}}}{k_{\mathrm{Pl}}}\right),$$
(32)
where $`F[k_{\mathrm{ph}}/k_{\mathrm{Pl}}]k_{\mathrm{ph}}/k_{\mathrm{Pl}}`$ in the large wavelength limit. Substituting the above form in (24), we get

$$I(\eta _{\mathrm{Pl}})=\frac{k_{\mathrm{Pl}}}{H}\times Q$$
(33)
where

$$Q=_{H\eta _i}^{H\eta _{\mathrm{Pl}}}\frac{d(H\eta )}{(H\eta )}F\left(\frac{k(H\eta )}{k_{\mathrm{Pl}}}\right)$$
(34)
Using the typical energy scale at the start of inflation to be GUT scale \[$`E10^{15}GeV`$, $`\eta _i10^{36}sec`$, $`H10^{11}GeV^1`$\], the prefactor $`k_{\mathrm{Pl}}/H`$ in the expression (33) is of the order $`10^{10}`$. This implies that unless the integral $`Q`$ is a small number (of the order $`10^9`$), the argument of the $`cosine`$ term in (29) rapidly oscillates and can be set to zero. \[This is because, in the finding $`C_l`$’s (spherical harmonics of the temperature anisotropy) we need to integrate over $`k`$ and the $`cosine`$ term vanishes.\]

Assuming the amount of inflation to be $`80e`$ foldings, we get $`\eta _{\mathrm{Pl}}10^{38}sec`$. To obtain a rough estimate of the integral $`Q`$, we consider four dispersions relations (II)– Padmanabhan’s, $`\kappa `$-Poincare, Unruh’s and Ahluwalia’s. The numerical value of the integral $`Q`$ in all these cases ranges between $`0.5`$ to $`0.7`$. \[For numerical work, we use MATHEMATICA owing to its convenience and flexibility it offers.\] Since, the argument of $`cosine`$ in the expression (29) is large (of the order $`10^9`$) for all practical purposes, it can be set to zero. Thus, Eq. (29) reduces to

$$k^3|C_k^{\mathrm{mod}}|^2=\frac{H^2}{2}+H^2|B_k|^2.$$
(35)
The following features are noteworthy regarding this result:

1. If the field modes are expanded as WKB modes and the initial state of the field is $`A_k=1;B_k=0`$, then the spectrum of fluctuations are flat scale-invariant for any general dispersion relation. \[In the case of Corley/Jacobson dispersion relation, $`\mathrm{\Omega }`$ becomes imaginary for $`k>k_{\mathrm{Pl}}`$ and hence the expansion of modes as WKB modes is not valid.\] In the standard inflationary scenario, the fluctuations are flat if the initial state of the field is $`a_k=1;b_k=0`$.
2. For an arbitrary initial state, the spectrum of fluctuations in the dispersive field theory models have a scale-invariant term plus a non-scale invariant term. The non-scale invariant term strongly depends on the initial state ($`B_k`$). However, the power spectrum does not depend strongly on the dispersion relation. Thus, the form of the dispersion relation does not play significant role in obtaining the corrections to the scale-invariant spectrum.
3. In our discussion, we have assumed $`|b_k|`$ and $`|B_k|`$ to be arbitrary functions of $`k`$. The form of $`|b_k|`$ and $`|B_k|`$ contribute to the non scale-invariant part of the power spectrum. Setting,

$$|B_k|^2=|b_k|^2+|b_k|(|b_k|^2+1)^{1/2},$$
(36)
we observe that the spectrum of fluctuations in the standard inflationary scenario and in dispersive field theory model are same. Hence, the signatures of the deviations from the flat scale-invariant spectrum from the CMBR observations due to quantum gravitational effects cannot be differentiated from the standard inflationary scenario with an arbitrary initial state. In other words, the signatures of the deviation from the flat scale-invariant spectrum is not necessarily due to the trans-Planckian effects. Even in the standard inflation, the deviations from the flat spectrum can be obtained for a different choice of vacuum state – other than the Bunch-Davies vacuum.

## III Conclusions

We have studied the quantum gravitational effects on the spectrum of primordial density perturbation in the inflationary epoch. We have incorporated the quantum gravitational corrections using non-linear dispersion relations to the quantum scalar field i. e., the scalar field satisfying the linear wave equation with higher spatial derivative terms propagating in the de Sitter space-time \[and ignored the back-reaction of the particles on the metric\]. We keep the form of the dispersion relation to be arbitrary and have not assumed any specific dispersive field theory model.

We have performed a general analysis of the power spectrum for a class of dispersion relations whose field modes can be expanded as WKB modes over a range of initial conditions for the quantum field. We have shown that the spectrum of fluctuations due to the quantum gravitational corrections has a scale-invariant term plus a non scale-invariant term which strongly depends on $`B_k`$ (the initial state) and the form of the dispersion relation $`\mathrm{\Omega }`$. However, the strong dependence on the dispersion relation cancels out due to the rapidly varying phase term. We have thus shown that the form of the dispersion relation do not play significant role in obtaining the corrections to the scale invariant spectrum. We have also shown that the deviations from the flat spectrum due to quantum gravitational effects (29) cannot be differentiated as opposed to the deviations from the flat spectrum for an arbitrary initial state in the standard inflationary scenario (13).

###### Acknowledgements.

The author would wish to thank T. Padmanabhan and K. Subramaniam for useful discussions. The author is being supported by the University Grants Commission, India through Senior Research Fellowship.


format: text

text:
# Charge fluctuations in a quark-antiquark system

## Abstract

It is pointed out that the recent data on charge fluctuations observed in heavy ion collisions are compatible with production of a system of weakly correlated constituent quarks and antiquarks.

It was recently suggested  (see also ) that measurements of the quantity

$$D=4\frac{<\delta Q^2>}{<N_{ch}>}ł1$$
(1)
can be used to distinquish the hadron gas in equilibrium from the quark-gluon plasma.

It was then argued in  that -after appropriate corrections for resonance production are taken into account- $`D3`$ for the hadron gas whereas for the quark-gluon plasma $`D1`$.

The preliminary data from CERES , NA49  and STAR  experiments indicate that the measured value of $`D`$ is close to that predicted for hadron gas and differs markedly from that expected for QGP, i.e. for a weakly correlated quark-gluon system.

This rises a question: what kind of structure of the partonic system produced at the early stage of collisions may be compatible with these experimental results? The present note is an attempt to answer this question.

The first observation is that the estimate of the numerator in (LABEL:1) is rather straigthforward. The argument (essentially the repetition of the argument given in ) goes as follows. Let us consider a system of several particle species (labelled by $`i`$) with charges $`q_i`$ and multiplicities $`n_i`$. Since

$`Q={\displaystyle \underset{i}{}}q_in_i<Q>={\displaystyle \underset{i}{}}q_i<n_i>ł2`$ (2)
we obtain

$`<\delta Q^2><Q^2><Q>^2={\displaystyle \underset{i}{}}(q_i)^2<n_i>+{\displaystyle \underset{i,k}{}}c_{ik}^{(2)}<n_i><n_k>q_iq_kł3`$ (3)
where $`c_{ik}^{(2)}`$ are the normalized two-particle correlation functions:

$`c_{ii}^{(2)}={\displaystyle \frac{<n_i(n_i1)>}{<n_i>^2}}1;c_{ik}^{(2)}={\displaystyle \frac{<n_in_k>}{<n_i><n_k>}}1ifik`$
$`ł4`$ (4)
If particles are uncorrelated the second term in (LABEL:3) vanishes. Moreover, one notes that the first term is the sum of only positive quantities, whereas in the second term cancellations are possible. For example, if all normalized two-particle correlations are identical, one obtains

$`<\delta Q^2>={\displaystyle \underset{i}{}}(q_i)^2<n_i>+c^{(2)}<Q>^2.ł5`$ (5)
and thus the second term disappers when $`<Q>=0`$.

From now on we shall thus restrict the discussion to the case when the second term can be neglected and thus

$`<\delta Q^2>={\displaystyle \underset{i}{}}(q_i)^2<n_i>ł6`$ (6)
For pion gas this means

$`<\delta Q^2>=<n_+>+<n_{}>=<N_{ch}>ł9`$ (7)
and thus

$`D=4.ł12`$ (8)
For a quark-gluon system we obtain

$`<\delta Q^2>={\displaystyle \frac{1}{9}}(4<n_u>+<n_d>+4<n_{\overline{u}}>+<n_{\overline{d}}>)ł10`$ (9)
(gluons of course do not contribute). In the simplest case when abundances of all quarks are identical we have

$`<\delta Q^2>={\displaystyle \frac{5}{18}}<N_q>ł11`$ (10)
where $`<N_q>`$ is the total number of quarks and antiquarks.

To calculate $`D`$ it is, however, also necessary to estimate $`<N_{ch}>`$ and this is of course model-dependent. The question is the role of gluons. If they decay into $`q\overline{q}`$ pairs before making hadrons, the number of hadrons may be very large. This is what happens in the QGP scenario. And this is the reason of the small value of $`D`$ obtained in . The argument is based on entropy. Qualitatively: since gluons have a large entropy, they have to produce many particles. Some of them will be charged and thus contribute to the denominator of (LABEL:1).

If one considers the $`q\overline{q}`$ coalescence scenario , however, the result is very different. In this scenario gluons are attached to quarks (forming constituent quarks) and thus do not provide important contribution to entropy. Hadrons are created simply by coalescence of the (constituent) $`q\overline{q}`$ pairs. Consequently, the total number of hadrons is about $`1/2`$ of the total number of quarks and antiquarks:

$`<N_h>={\displaystyle \frac{1}{2}}<N_q>ł13`$ (11)
It is natural to assume that $`2/3`$ of them are charged, so that we obtain

$`<N_{ch}>={\displaystyle \frac{2}{3}}<N_h>={\displaystyle \frac{1}{3}}<N_q>.ł14`$ (12)
and thus

$`D={\displaystyle \frac{10}{3}}=3.333ł15`$ (13)
which is only slightly smaller than (LABEL:12).

The values of $`D`$ given by (LABEL:12) and (LABEL:15) should still be corrected for resonance production . Since the difference between them is fairly small, however, one should not expect a great difference between the corrected results, either.

We thus conclude that the data on charge fluctuations are not incompatible with the idea that a (quasi)uncorrelated constituent quark system is created before hadronization. This observation strenghtens earlier arguments leading to the coalescence picture of hadronization (see e.g.  -).

Two comments are in order.

(i) The present data are given for fairly small rapidity intervals and thus important effects of charge transport between the boundaries of the considered phase-space region may be present , destroying the argument of Refs.  and . This effect may perhaps also be responsible for the rather large values of $`D`$ (exceeding 4) observed at SPS energies . Data for larger intervals would be very useful to sort out this question.

(ii) An argument similar to that presented here can be also applied to fluctuations of the net baryon number . The formulae (LABEL:2) - (LABEL:6) are valid with the substitution of the charges by baryon numbers. In particular, for the hadron gas we have

$`<\delta B^2>={\displaystyle \underset{i}{}}B_i^2<n_i>=<N_B>ł16`$ (14)
whereas for a parton system one obtains

$`<\delta B^2>={\displaystyle \frac{1}{9}}<N_q>.ł17`$ (15)
Thus the ratio $`<\delta B^2>/<N_{ch}>0.1`$ for the QGP while it is expected to be about three times larger in the coalescence model. Again, the corrections for resonance decays must be included to obtain a realistic estimate.

(iii) The constituent quark coalescence model discussed here can also be intepreted as a quark-gluon plasma far from equilibrium. This poses an apparent problem, because it has been shown in numerous papers  that the equilibrium scenario agrees very well with the data, in particular with those on particle ratios (chemical equilibrium). As I have already mentioned, however, particle ratios were analysed in detail by the Budapest group  and the coalescence model was also found to be in excellent agreement with data. I thus feel that it deserves a serious attention. Clearly, however, more work is needed to clarify fully the situation.

Acknowledgements

I would like to thank K.Fialkowski, M.Jezabek, V.Koch and M.Nowak for discussions. This investigation was supported in part by Subsydium of Foundation for Polish Science NP 1/99.


text:
# Ground State Vortex Lattice Structures in 𝑑-wave Superconductors

## I Introduction

An external magnetic field enters a (type II) superconductor as a collection of quantized magnetic flux tubes. The flux tubes with associated supercurrent vortices form a triangular lattice  as is seen in conventional superconductors . Deviations from this structure are of considerable interest. In conventional superconductors, the observed deviations have been attributed to the anisotropy of the underlying one electron energy spectrum. In heavy fermion and high Tc superconductors, an additional and potentially very interesting reason is the existence of an unconventional superconducting order parameter, with gaps which have nodes and change sign. Indeed there is experimental evidence both in heavy fermion systems  the exotic superconductor SrRuO4  and in cuprate superconductors that the vortex lattice is not triangular (for some field and temperature regimes). Recently, measurements of small angle neutron diffraction from untwinned YBa2Cu3O7-δ single crystals shows a very well formed triangular lattice that undergoes an orientational transition from along $`a`$ axis to along $`b`$ axis for a $`3T`$ magnetic field at $`33^o`$ to the $`c`$-axis. The reasons for possible nontriangular structure as well as for the structural transition are not fully established, and are the subject of considerable theoretical work . In high Tc superconductors, entropic effects, abetted by high transition temperature as well as weak interlayer coupling play an important role, and the classical statistical mechanics of interacting, meandering flux lines, of the vortex fluid phase, and the solid fluid transition has developed into a major theoretical and experimental subfield .

In this paper, we focus on the nature of the flux lattice at temperatures well below the superconducting transition, where vortex configurational entropy effects mentioned above are negligible. The ground state structures and structural transition then directly reflect the electronic peculiarities of the superconductor, and thus probe the latter. For cuprate superconductors, a number of measurements show that the superconducting gap $`\mathrm{\Delta }_𝐤`$ has nodes , has a magnitude with a $`(\mathrm{\Delta }_0/2)|\mathrm{cos}k_xa\mathrm{cos}k_ya|`$ dependence on the two dimensional wave vector $`𝐤`$ across the Fermi surface  and that transport lifetimes of quasiparticles are long  for $`TT_c`$. Thus one can assume well defined low energy, nodal quasiparticles, with an experimentally determined one electron dispersion $`ϵ_𝐤`$  and gap function $`\mathrm{\Delta }_𝐤`$ . The question of interest is the effect of the zero gap, anisotropic Dirac like linear quasiparticle excitation spectrum on the interaction between vortices, and thus on vortex lattice structure. There is considerable evidence e.g., from magnetic field dependent electronic specific heat , electronic thermal conductivity , and superfluid density  measurements that an external magnetic field going in as vortices has a strong effect on electronic states, changing their density and lifetime. The relevant issue here is somewhat the reverse, namely the effect of the quasiparticles on interaction between vortices. The order parameter phase associated with the vortex, and the related magnetic vector potential together constitute the superfluid velocity field $`𝐯_s(𝐫)[=_l𝐯_s(𝐫𝐑_l)]`$ where the vortices are located at points $`𝐑_l`$. The extra superfluid kinetic energy, being quadratic in $`v_s`$, clearly has a part that depends on two vortex coordinates and is thus structure sensitive. In addition to this ‘diamagnetic’ term, which is the additional kinetic energy of the rigidly moving superfluid and which is minimized for a triangular lattice , there is another ‘paramagnetic’ term due to the perturbation of quasiparticles by the superfluid current via the term $`𝐤𝐯_s`$. This causes virtual particle hole excitations; two vortices interact via the exchange of quasiparticle quasihole pairs. This polarization term depends on the quasiparticle excitation spectrum. In a clean $`s`$-wave superconductor, the process leads to an additional isotropic interaction between vortices of order $`(H/H_{c_2})`$ relative to the diamagnetic term. However, in a $`d`$-wave superconductor where the excitation gap vanishes along some (nodal) directions, one expects the nonlocal polarizeability to be larger as well as anisotropic; this gives rise to an interaction between vortices which depends on the orientation of the line joining them with the respect to crystaline axes, and consequently can be the cause of novel long range positional order.

The ground state energy arising from quasiparticle / hole mediated interaction between vortices depends linearly on the nonlocal current susceptibility $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$, for wave vectors $`𝐪`$ equal to the reciprocal lattice vectors $`𝐆`$ of the vortex lattice. Because the gap as well as density of quasiparticle states vanish linearly near the node, $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$ is proportional to $`|q_x|`$ or $`|q_y|`$ for small $`q`$. This nonanalytic behaviour, noticed first by Kosztin and Leggett , has also been discussed by Franz et al  who were the first to analyze microscopically its effect, as well as of the anisotropy in $`\chi _{\alpha \beta }^\mathrm{p}`$, on vortex lattice structure at $`T=0`$. These authors found a rich phase diagram in the field temperature plane, with a centered rectangular lattice at $`T=0`$ whose inner angle varies continuously as a function of field, as well as a sudden orientational transition at higher temperatures, and a transition to a centered square lattice for very high fields and low $`T`$. In obtaining these results, Franz et al made a “local” approximation for the gap function, i.e., assumed $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_{𝐤+𝐆}`$, and more importantly they assumed a momentum independent quasiparticle current which leads to a response function $`\chi _{xx}^\mathrm{p}(𝐆)=\chi _{yy}^\mathrm{p}(𝐆)`$. The anisotropy then enters only through $`\chi _{xy}^\mathrm{p}(𝐆)`$. We carry out here a more detailed and realistic calculation of the nonlocal susceptibility, considering the strong $`𝐤`$ dependence of the gap function $`\mathrm{\Delta }_𝐤`$ and quasiparticle current $`𝐣_𝐤`$ properly, and using a realistic one electron dispersion. The diagonal terms $`\chi _{xx}^\mathrm{p}`$ and $`\chi _{yy}^\mathrm{p}`$ are unequal and large, and this anisotropy is seen to be the underlying cause of the transition. The contibution of the off diagonal susceptitiblity $`\chi _{xy}^\mathrm{p}`$ is smaller than that of the diagonal susceptibilities. Our results for structural stability (at $`T=0`$) are therefore quite different from those of Franz et al .

Confining ourselves to $`T=0`$, we find, as summarized in a phase diagram (see Fig. 1), that the stable lattice at low fields is triangular. At about 5 Tesla (for the parameters chosen) the orientation of the smallest $`𝐆`$ vector changes from being along one of the axes to being along the order parameter node direction, because the system is most susceptible to excitations with wave vector along the node. We have analyzed the driving force for this transition, both analytically and numerically, and find that it arises from a subtle balance between the term linear in $`|𝐆|`$ and the quadratic term, which are slightly different for the two orientations. The field scale for the transition is approximately given by the condition $`(taG/\mathrm{\Delta }_0)1`$ which is natural on dimensional grounds. The Fermi velocity is $`ta`$ and the energy scale associated with the superflow $`(\theta )`$ with Fourier component $`G`$ is therefore $`taG`$. The polarizability or susceptibility has an energy scale $`(1/\mathrm{\Delta }_0)`$ where $`\mathrm{\Delta }_0`$ is the gap magnitude which sets the scale for quasiparticle excitation energies. Thus the dimensionless susceptibility of interest is $`(taG/\mathrm{\Delta }_0)`$. For realistic parameters $`t,a`$, and $`\mathrm{\Delta }_0`$ this translates $`[(taG_c/\mathrm{\Delta }_0)0.37]`$ to a field scale of 5.2T.

We find that the node oriented triangular lattice is stable till about 28T, whereupon a discontinuous transition to a centered square lattice takes place. This structure, which is orientationally commensurate with the symmetry of the quasiparticle dispersion, is probably the most stable $`T=0`$ phase when electronic commensurate effects dominate. However, the calculated field scale is large enough that the London approximation used, valid for $`HH_{c_2}`$, is not reliable, because vortex core effects cannot be neglected at these high fields.

In the next section (Section II), we describe the model and the theoretical approach used. The tight binding quasiparticle Hamiltonian is decomposed into an unperturbed part $`H_0`$ and a term $`H_I`$ due to the quasiparticle vortex interaction. The free energy or the ground state energy can be obtained as a power series in $`H_I`$ or equivalently the density of vortices. For low vortex densities $`(HH_{c_2})`$ the leading or $`n_v^2`$ term is sufficient, and describes quasiparticle-hole mediated vortex interaction, in addition to the bare superfluid kinetic energy. We discuss the former carefully in terms of the nonlocal, anisotropic current susceptibility $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$ since the energy can be expressed as a reciprocal lattice vector sum over $`\chi _{\alpha \beta }^\mathrm{p}(𝐆)`$. We obtain $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$ semianalytically for small $`q`$ at $`T=0`$, as well as numerically (Section III).

The calculations for different two dimensional structures are discussed in section III. For a given magnetic field $`B`$ the most general centered rectangular lattice \[$`a_1,a_2`$\] can be described in terms of an angle $`\theta `$ related to the aspect ratio $`(a_1/a_2)`$ as $`\mathrm{tan}\theta =a_1/a_2`$, and an orientation $`\varphi `$ with respect to crystal axes (Fig. 2). We compute the ground state energy as a function of these two variables for different magnetic fields. The basic vortex related electronic energy parameters are the following. The vortex has three energy scales, namely the diamagnetic single vortex energy, of order 3450$`K`$ per vortex, the vortex interaction energy, the diamagnetic part of which has a value $`1440K`$ (for nearest neighbor vortices), and the paramagnetic interaction term which is about $`350K`$ at 5T field. The last, and the smallest term is structure sensitive and is of interest here. The ground state is analyzed as a function of $`\theta `$, $`\varphi `$ for several field values in section III. It turns out that the structure sensitive part of the last (paramagnetic vortex interaction) term is extremely small, of order a few degrees per vortex. This has obvious implications for the observability of the transition, because the structural changes predicted and the clean limit anisotropies obtained can be easily overwhelmed by effects of disorder, eg., vortex pinning and the muting of the paramagnetic susceptibility anisotropy, and nonanalyticity by disorder. However, the size of the structure sensitive terms is larger, the greater the $`(v_F/v_\mathrm{\Delta })`$ ratio or anisotropy. One can thus imagine situations where this effect is quite large.

In the concluding section (Section IV) we briefly discuss thermal effects, the consequence of the predicted transition and their observability, the reason for which our result doffers from earlier, and the experimentally observed structural transitions.

## II Theory

### A Model

We consider a two dimensional lattice model with nearest neighbor and next nearest neighbor hopping for a CuO2 plane of high $`T_c`$ superconductors. The (mean field) Hamiltonian in this model is given by

$$H_0=t\underset{<ij>_1\sigma }{}(c_{i\sigma }^{}c_{j\sigma }+h.c.)+t^{}\underset{<<ij>>_1\sigma }{}(c_{i\sigma }^{}c_{j\sigma }+h.c.)+\underset{<ij>}{}(\mathrm{\Delta }_{ij}c_i^{}c_j^{}+h.c.)\mu \underset{i}{}c_{i\sigma }^{}c_{i\sigma },$$
(1)
where $`t`$ and $`t^{}`$ are nearest neighbor and next nearest neighbor hopping integrals respectively. This corresponds for appropriate choices of $`t`$ and $`t^{}`$ to an open Fermi surface which is observed in angle resolved photo emission experiments. The pair amplitude $`\mathrm{\Delta }_{ij}`$ is considered to be $`d_{x^2y^2}`$-wave like, i.e., $`\mathrm{\Delta }_{i,i\pm a\widehat{x}}=\mathrm{\Delta }_{i,i\pm a\widehat{y}}`$, where $`a`$ is the lattice constant in a square lattice, and $`\mu `$ is the chemical potential.

When we apply a magnetic field beyond the lower critical field $`H_{c_1}`$ in high $`T_c`$ superconductors, the magnetic field goes into the system in the form of vortices. The magnetic induction is screened over a length $`\lambda `$, the penetration depth. The pair amplitude acquires a phase: $`\mathrm{\Delta }_{ij}\mathrm{\Delta }_{ij}\mathrm{exp}[i\theta _{ij}]`$, where $`\theta _{ij}`$ is the sum of polar angles of all the vortices measured with respect to a particular axis, for the centre of mass of the pair $`ij`$. We write $`\theta _{ij}`$ as $`(\theta _i+\theta _j)/2`$, (as an average of the angles of individual coordinates of the Cooper pairs), which is consistent upto $`𝒪(1/k_F\xi )^2`$, where $`k_F`$ is the Fermi momentum and $`\xi `$ is the superconductive coherence length. The pair amplitude (order parameter magnitude) vanishes at the center of the core of a vortex, and over a distance $`\xi `$ it acquires its uniform value. For a collection of vortices with $`HH_{c_2}`$, i.e., with intervortex spacing $`\xi `$, or the London limit, we assume the order parameter magnitude be uniform, throughout the superconductor (there are $`\delta `$ function sources of phase rotation at the locations of vortices). There is a vector potential $`𝐀`$ such that $`\times 𝐀(𝐫)=B(𝐫)`$ where $`B(𝐫)`$ is the local magnetic induction, along the $`c`$-axis. Its effect in this model is to change the hopping integrals to

$$(t,t^{})(t,t^{})\mathrm{exp}[i(e/\mathrm{}c)_{𝐫_i}^{𝐫_j}𝐀𝑑𝐥]$$
(2)
for hopping from site $`j`$ to site $`i`$. We then make a gauge transformation: $`c_{i\sigma }c_{i\sigma }e^{i\theta _i/2}`$. We thus obtain the Hamiltonian as

$`H`$ $`=`$ $`t{\displaystyle \underset{<ij>_1\sigma }{}}(c_{i\sigma }^{}c_{j\sigma }\mathrm{exp}[i(\theta _i\theta _j)/2+i(e/\mathrm{}c){\displaystyle _{𝐫_i}^{𝐫_j}}𝐀d𝐥]+h.c.)+{\displaystyle \underset{<ij>_1\sigma }{}}(\mathrm{\Delta }_{ij}c_i^{}c_j^{}+h.c.)`$ (4)
$`+t^{}{\displaystyle \underset{<<ij>>}{}}(c_{i\sigma }^{}c_{j\sigma }\mathrm{exp}[i(\theta _i\theta _j)/2+i(e/\mathrm{}c){\displaystyle _{𝐫_i}^{𝐫_j}}𝐀d𝐥]+h.c.)\mu {\displaystyle \underset{i}{}}c_{i\sigma }^{}c_{i\sigma }.`$
The phase difference between two nearest or next nearest neighbour sites can be expressed as

$$\frac{1}{2}(\theta _i\theta _j)+(e/\mathrm{}c)_{𝐫_i}^{𝐫_j}𝐀𝑑𝐥(𝐫_i𝐫_j)(\frac{1}{2}_i\theta (e/\mathrm{}c)𝐀_i)(m/\mathrm{})(𝐫_i𝐫_j)𝐯_s(𝐫_i).$$
(5)
Here the superfluid velocity

$$𝐯_s(𝐫)=\frac{1}{m}\left[\frac{\mathrm{}}{2}\theta \frac{e}{c}𝐀(𝐫)\right]$$
(6)
for a single vortex and for a collection of vortices, $`𝐯_s(𝐫)=_l𝐯_s(𝐫𝐑_l)`$ (where the vortices are located at $`𝐑_l`$. We then assume that the phase difference between two neighboring lattice sites is very small (which is certainly true in the London limit) so that we expand exponentials in Eq. (4) upto quadratic terms.

Using Eqs. (5) and (4) in Eq.(3) for $`H`$, and expanding upto quadratic order in the small quantity $`𝐯_s(𝐫)`$ we have

$$H=H_0+H_I+H_{II},$$
(7)
where the free Hamiltonian

$$H_0=\underset{𝐤,\sigma }{}\xi _𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }+\underset{𝐤}{}[\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+h.c.],$$
(8)
with $`\xi _𝐤=2t[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)]+4t^{}\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)\mu `$ and $`\mathrm{\Delta }_𝐤=(\mathrm{\Delta }_0/2)[\mathrm{cos}(k_xa)\mathrm{cos}(k_ya)]`$, $`\mathrm{\Delta }_0`$ being the maximum quasiparticle excitation gap. Here $`𝐤`$ lies in the first atomic Briliouin zone, i.e., $`\frac{\pi }{a}(k_x,k_y)\frac{\pi }{a}`$. A typical structure of the Fermi surface is shown in Fig. 3. Gapless quasiparticle excitations exist along $`k_x=\pm k_y`$ directions as noted in the figure. The interaction term (first order in $`v_s`$) can now be expressed as

$`H_I`$ $`=`$ $`2(at/\mathrm{}){\displaystyle \underset{𝐤}{}}{\displaystyle \underset{𝐆>0}{}}[c_{𝐤\sigma }^{}c_{𝐤+𝐆\sigma }c_{𝐤+𝐆\sigma }^{}c_{𝐤\sigma }][mv_s^x(𝐆)\mathrm{sin}(k_xa)+mv_s^y(𝐆)\mathrm{sin}(k_ya)]`$ (11)
$`4(at^{}/\mathrm{}){\displaystyle \underset{𝐤}{}}{\displaystyle \underset{𝐆>0}{}}[c_{𝐤\sigma }^{}c_{𝐤+𝐆\sigma }c_{𝐤+𝐆\sigma }^{}c_{𝐤\sigma }][mv_s^x(𝐆)\mathrm{sin}(k_xa)\mathrm{cos}(k_ya)`$
$`+mv_s^y(𝐆)\mathrm{sin}(k_ya)\mathrm{cos}(k_xa)]`$
$``$ $`{\displaystyle \underset{𝐤\sigma }{}}{\displaystyle \underset{G>0}{}}V_{𝐤,𝐆}[c_{𝐤\sigma }^{}c_{𝐤+𝐆\sigma }c_{𝐤+𝐆\sigma }^{}c_{𝐤\sigma }]`$ (12)
Here $`V_{𝐤,𝐆}`$ is purely imaginary. The term $`H_{II}`$ is quadratic in $`v_s`$ and contributes to the free energy as a diamagnetic term. It is given by

$$H_{II}=2N_s^0(tt^{})(a^2/\mathrm{}^2)\underset{𝐆}{}mv_s^\alpha (𝐆)mv_s^\alpha (𝐆)$$
(13)
with $`N_s^0`$ being the number of superfluid carriers, and $`\alpha `$ refers to cartesian variables $`x`$ and $`y`$. Paramekanti et al.  have recently shown that the quantum phase fluctuation of the order parameter reduces the superfluid density considerably. We thus reexpress the term $`H_{II}`$ phenomenologically in terms of the measured $`\lambda `$ as

$$H_{II}=\frac{dA}{2\lambda ^2}\left(\frac{c^2}{4\pi e^2}\right)\underset{𝐆}{}mv_s^\alpha (𝐆)mv_s^\alpha (𝐆),$$
(14)
where $`d`$ is the mean interlayer separation of weakly coupled superconducting layers and $`A`$ is the area of the system.

### B Free Energy

We now calculate the free energy as a power series in $`𝐯_s(𝐫)`$ or equivalently the vortex density. The diamagnetic or Ginzburg-Landau term, of first order in $`H_{II}`$, is the largest, and the structure sensitive part of it is known to be minimized for a triangular lattice (1). The energy does not depend on its orientation with respect to the crystal lattice. We are interested here additionally in the paramagnetic term, of second order in $`H_I`$. Including this, and the magnetic field energy contribution, the free energy to second order in vortex density is given (per unit length along the $`c`$ axis) by

$$\mathrm{\Delta }\mathrm{\Omega }=\frac{1}{2Ad\mathrm{}^2}\underset{𝐆}{}mv_s^\alpha (𝐆)\left[\chi ^d\delta _{\alpha \beta }\chi _{\alpha \beta }^\mathrm{p}(𝐆)\right]mv_s^\beta (𝐆)+\frac{1}{8\pi A}\underset{𝐆}{}B_𝐆B_𝐆,$$
(15)
where $`𝐆`$ is the reciprocal vector of the vortex lattice. The individual vortex energy is not included here, as it is not relevant for the question of vortex lattice structure. $`\chi ^d=(c^2\mathrm{}^2d/4\pi e^2\lambda ^2)`$ is the diamagnetic term arising from the term $`H_{II}`$ (Eq.14) to first order, and $`\chi _{\alpha \beta }^\mathrm{p}(𝐆)`$ is the paramagnetic current susceptibility due to second order contribution from $`H_I`$. Higher order contributions are neglected since the expansion parameter is $`(n_v/n)`$ where $`n_v`$ is the vortex density and $`n`$ is the electron density. This ratio is obviously much smaller than one.

The paramagnetic susceptibility $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$ is expressed as

$$\chi _{\alpha \beta }^\mathrm{p}(𝐪)=\frac{1}{(2\pi )^2}_{\pi /a}^{\pi /a}𝑑k_x_{\pi /a}^{\pi /a}𝑑k_y\gamma _{\alpha \beta }(𝐤)\mathrm{\Pi }(𝐤,𝐪).$$
(16)
The current operator dependent terms $`\gamma _{\alpha \beta }(𝐤)`$ are explicitly given as

$`\gamma _{xx}(𝐤)`$ $`=`$ $`\left[2a\mathrm{sin}(k_xa)(t2t^{}\mathrm{cos}(k_ya))\right]^2,`$ (18)
$`\gamma _{yy}(𝐤)`$ $`=`$ $`\left[2a\mathrm{sin}(k_ya)(t2t^{}\mathrm{cos}(k_xa))\right]^2,`$ (19)
$`\gamma _{xy}(𝐤)`$ $`=`$ $`\left[2a\mathrm{sin}(k_xa)(t2t^{}\mathrm{cos}(k_ya))\right]\left[2a\mathrm{sin}(k_ya)(t2t^{}\mathrm{cos}(k_xa))\right],`$ (20)
$`=`$ $`\gamma _{yx}(𝐤).`$ (21)
The zero frequency susceptibility of wave vector $`𝐪`$ for quasiparticle quasihole of momentum $`𝐤`$ is $`\mathrm{\Pi }(𝐤,𝐪)`$ and has the form

$$\mathrm{\Pi }(𝐤,𝐪)=\frac{1}{E_𝐤+E_{𝐤+𝐪}}\left[1\frac{\xi _𝐤\xi _{𝐤+𝐪}+\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}}{E_𝐤E_{𝐤+𝐪}}\right],$$
(22)
with the quasiparticle energy $`E_𝐤=\sqrt{\xi _𝐤^2+\mathrm{\Delta }_𝐤^2}`$.

It is expected that the above susceptibilities are anisotropic due to the nonlocal nature of $`\mathrm{\Delta }_{ij}`$, reflected in the $`𝐤`$ dependence of $`\mathrm{\Delta }_𝐤`$. Though anisotropic, they possess certain symmetries: $`\chi _{\alpha \alpha }^\mathrm{p}(q_x,q_y)=\chi _{\alpha \alpha }^\mathrm{p}(q_x,q_y)=\chi _{\alpha \alpha }^\mathrm{p}(q_x,q_y)`$; $`\chi _{xy}^\mathrm{p}(q_x,q_y)=\chi _{xy}^\mathrm{p}(q_x,q_y)=\chi _{xy}^\mathrm{p}(q_x,q_y)`$, and $`\chi _{xx}^\mathrm{p}(q_x,q_y)=\chi _{yy}^\mathrm{p}(q_y,q_x)`$. These symmetries suggest that the susceptibilities are functions of $`|q_x|`$, $`|q_y|`$ and sign($`q_xq_y`$) only. A naive perturbative expansion of $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$ in powers of $`q`$ fails since the coefficient of quadratic term in $`q`$ is divergent, due to the vanishing of $`\mathrm{\Delta }_𝐤`$ on the Fermi surface at $`k_x=\pm k_y`$ points. We however proceed to evaluate these analytically as follows.

We write

$$\chi _{\alpha \beta }^\mathrm{p}(q_x,q_y)=\underset{j=1}{\overset{4}{}}\chi _{\alpha \beta }^{\mathrm{p},j}(q_x,q_y),$$
(23)
where $`\chi _{\alpha \beta }^{\mathrm{p},j}`$ is the contripution of $`j`$-th quadrant $`(j=`$1–4) of $`k`$-space. For instance,

$$\chi _{\alpha \beta }^{\mathrm{p},1}(q_x,q_y)=\frac{1}{(2\pi )^2}_0^{\pi /a}𝑑k_x_0^{\pi /a}𝑑k_y\gamma _{\alpha \beta }(𝐤)\mathrm{\Pi }(𝐤,𝐪)$$
(24)
is the contribution due to 1st quadrant. We present the calculation of $`\chi _{\alpha \beta }^{\mathrm{p},1}(q_x,q_y)`$ below in detail.

In terms of an alternative coordinate system $`(k_1,k_2)`$ whose origin is at the nodal point on the Fermi surface as shown in Fig. 4, the old coordinates in the first quadrant are expressed as $`k_x=\frac{1}{\sqrt{2}}(k_0+k_1k_2)`$ and $`k_y=\frac{1}{\sqrt{2}}(k_0+k_1+k_2)`$, where $`k_0`$ is defined as $`\mu =4t\mathrm{cos}(k_0a/\sqrt{2})+4t^{}\mathrm{cos}^2(k_0a/\sqrt{2})`$. We use $`\xi _𝐤\mathrm{}v_Fk_1`$ and $`\mathrm{\Delta }_𝐤\mathrm{}v_\mathrm{\Delta }k_2`$ in linear form, where the Fermi velocity $`v_F=(\frac{4a}{\mathrm{}\sqrt{2}})[t2t^{}\mathrm{cos}\frac{k_0a}{\sqrt{2}}]\mathrm{sin}\frac{k_0a}{\sqrt{2}}`$ and $`v_\mathrm{\Delta }=(\frac{a}{\mathrm{}\sqrt{2}})\mathrm{\Delta }_0\mathrm{sin}\frac{k_0a}{\sqrt{2}}`$ is the velocity of quasiparticles along the $`k_2`$ direction. Since in $`d`$-wave superconductors $`v_Fv_\mathrm{\Delta }`$, the phase space of $`k_2`$ effectively is much larger than that of $`k_1`$ for a given value of quasiparticle energy. We observe that this is the cause of strong anisotropy in the diagonal susceptibilities as we see below. If $`\varphi `$ be the angle of a $`𝐤`$-vector with $`k_1`$ axis in this new coordinate system, $`\xi _𝐤E_𝐤\mathrm{cos}\varphi `$ and $`\mathrm{\Delta }_𝐤E_𝐤\mathrm{sin}\varphi `$. By Taylor expansion in $`q`$ we find $`\xi _𝐤\xi _{𝐤+𝐪}+\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}E_𝐤(E_𝐤+\alpha _𝐪)`$ and $`E_𝐤E_{𝐤+𝐪}|E_𝐤^2+E_𝐤\alpha _𝐪+\beta _𝐪^2|`$, where $`\alpha _𝐪=\mathrm{}(q_1v_F\mathrm{cos}\varphi +q_2v_\mathrm{\Delta }\mathrm{sin}\varphi )`$ and $`\beta _𝐪^2=(\mathrm{}^2/2)(q_1v_F\mathrm{sin}\varphi q_2v_\mathrm{\Delta }\mathrm{cos}\varphi )^2`$ with $`q_{1,2}=(q_y\pm q_x)/\sqrt{2}`$ respectively. Since $`\alpha _𝐪`$ is negative for some region of $`\varphi `$, the quantity $`(E_𝐤^2+E_𝐤\alpha _𝐪+\beta _𝐪^2)`$ may be negative as well as positive which we refer below as the region I and II respectively. It is negative in the regime $`E_1^0<E_𝐤<E_2^0`$, where $`E_1^0(\beta _𝐪^2/\alpha _𝐪)`$ and $`E_2^0(\alpha _𝐪+\beta _𝐪^2/\alpha _𝐪)`$. Expanding $`\gamma _{\alpha \beta }(𝐤)`$ upto linear in $`E_𝐤`$, we perform the integrals over $`E_𝐤`$ in for both regions I and II separately to obtain

$`\chi _{xx}^{\mathrm{p},1}(q_x,q_y){\displaystyle \frac{a^2}{\pi ^2\mathrm{}^2v_Fv_\mathrm{\Delta }}}`$ $`[{\displaystyle _I}d\varphi \{_1(\alpha _𝐪+2(\beta _𝐪^2/\alpha _𝐪)\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)`$ (28)
$`+{\displaystyle \frac{2}{3}}(𝒟_\varphi ^+_2+𝒟_\varphi ^{}_3)(\alpha _𝐪^23\beta _𝐪^23\beta _𝐪^2\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)\}`$
$`+{\displaystyle _{II}}d\varphi \{2\mathrm{ln}(2)(\beta _𝐪^2/\alpha _𝐪)_12\beta _𝐪^2(𝒟_\varphi ^+_2+𝒟_\varphi ^{}_3)\mathrm{ln}(\alpha _𝐪/\mathrm{\Delta }_0)\}],`$
$`\chi _{yy}^{\mathrm{p},1}(q_x,q_y){\displaystyle \frac{a^2}{\pi ^2\mathrm{}^2v_Fv_\mathrm{\Delta }}}`$ $`[{\displaystyle _I}d\varphi \{_1(\alpha _𝐪+2(\beta _𝐪^2/\alpha _𝐪)\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)`$ (31)
$`+{\displaystyle \frac{2}{3}}(𝒟_\varphi ^+_3+𝒟_\varphi ^{}_2)(\alpha _𝐪^23\beta _𝐪^23\beta _𝐪^2\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)\}`$
$`+{\displaystyle _{II}}d\varphi \{2\mathrm{ln}(2)(\beta _𝐪^2/\alpha _𝐪)_12\beta _𝐪^2(𝒟_\varphi ^+_3+𝒟_\varphi ^{}_2)\mathrm{ln}(\alpha _𝐪/\mathrm{\Delta }_0)\}],`$
$`\chi _{xy}^{\mathrm{p},1}(q_x,q_y){\displaystyle \frac{a^2}{\pi ^2\mathrm{}^2v_Fv_\mathrm{\Delta }}}`$ $`[{\displaystyle _I}d\varphi \{_1(\alpha _𝐪+2(\beta _𝐪^2/\alpha _𝐪)\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)`$ (34)
$`+{\displaystyle \frac{2}{3}}(_2+_3)\mathrm{cos}\varphi (\alpha _𝐪^23\beta _𝐪^23\beta _𝐪^2\mathrm{ln}|\alpha _𝐪/\mathrm{\Delta }_0|)\}`$
$`+{\displaystyle _{II}}d\varphi \{2\mathrm{ln}(2)(\beta _𝐪^2/\alpha _𝐪)_12\beta _𝐪^2(_2+_3)\mathrm{cos}\varphi \mathrm{ln}(\alpha _𝐪/\mathrm{\Delta }_0)\}],`$
where $`𝒟_\varphi ^\pm =\mathrm{cos}\varphi \pm (v_F/v_\mathrm{\Delta })\mathrm{sin}\varphi `$,

$`_1`$ $`=`$ $`(t^2+\mu t^{})\mathrm{sin}^2(k_0a/\sqrt{2}),`$ (36)
$`_2`$ $`=`$ $`2\sqrt{2}({\displaystyle \frac{a}{\mathrm{}v_F}})t^{}(t2t^{}\mathrm{cos}(k_0a/\sqrt{2}))\mathrm{sin}^3(k_0a/\sqrt{2}),`$ (37)
$`_3`$ $`=`$ $`\sqrt{2}({\displaystyle \frac{a}{\mathrm{}v_F}})(t^2+\mu t^{})\mathrm{sin}(k_0a/\sqrt{2})\mathrm{cos}(k_0a/\sqrt{2}),`$ (38)
and $`_I`$ and $`_{II}`$ represent the integrals over $`\varphi `$ $`(0\varphi 2\pi )`$ for the regime of $`\varphi `$ in which $`\alpha _𝐪<0`$ and $`>0`$ respectively. We, similarly, calculate $`\chi _{\alpha \beta }^{\mathrm{p},j}`$ for three other quadrants. We observe that $`\chi _{xx}^{\mathrm{p},1}`$ and $`\chi _{yy}^{\mathrm{p},1}`$ differ substantially through the last terms within both the integrals $`_I`$ and $`_{II}`$ in their expressions (28) and (31) since $`v_Fv_\mathrm{\Delta }`$. These lead to an anisotropic diagonal susceptibility. We note that these terms also the terms involving $`_{II}`$ arise due to keeping the linear dependences of $`k_1`$ and $`k_2`$ in $`\gamma _{\alpha \beta }(𝐤)`$. However for an approximation $`k_x=k_y=k_0/\sqrt{2}`$ in $`\gamma _{\alpha \beta }(𝐤)`$, $`\chi _{xx}^\mathrm{p}=\chi _{yy}^\mathrm{p}`$ as obtained by Franz et al .

Since the angular integrals in the expressions for $`\chi _{\alpha \beta }^\mathrm{p}`$ cannot be performed analytically, we numerically integrate these to obtain $`\chi _{\alpha \beta }^\mathrm{p}`$ in the next section. We shall then compare these semianalytically obtained $`\chi _{\alpha \beta }^\mathrm{p}`$ with those completely numerically obtained through Eqs.(16)–(22).

We now turn to obtain the equation for the vector potential in the gauge $`𝐆𝐀_𝐆=0`$ (from Eq. (15)). By minimizing the energy with respect to as $`𝐀_𝐆`$, we have

$$(A_𝐆)_\alpha =\frac{4\pi e}{\mathrm{}cG^2d}[\chi ^d\delta _{\alpha \beta }\chi _{\alpha \beta }^\mathrm{p}(𝐆)]\left[\frac{1}{2}(\theta )_𝐆\frac{e}{\mathrm{}c}A_𝐆\right]_\beta .$$
(39)
We thus obtain

$$\left[\frac{\mathrm{}}{2}(\theta )_𝐆\right]_\alpha =\left[G^2Q_{\beta \alpha }^1+\delta _{\beta \alpha }\right]\left(\frac{e}{c}A_𝐆\right)_\beta $$
(40)
where

$$Q_{\alpha \beta }(𝐪)=\frac{1}{\lambda ^2}\delta _{\alpha \beta }\left(\frac{4\pi e^2}{c^2\mathrm{}^2d}\right)\chi _{\alpha \beta }^\mathrm{p}(𝐪)$$
(41)
Using Eq. (40) in Eq. (15), we get

$$\mathrm{\Delta }\mathrm{\Omega }=\frac{1}{8\pi A}\underset{𝐆}{}B_𝐆\left[1+\frac{G_\alpha Q_{\alpha \beta }G_\beta }{\text{Det}Q}\right]B_𝐆.$$
(42)
We now express $`B_𝐆`$ in terms of $`N_v,\mathrm{\Phi }_0,Q_{\alpha \beta }`$ and $`𝐆`$. For a vortex lattice,

$$(\theta )_𝐆=2i\pi N_v\frac{𝐆\times \widehat{e}_z}{G^2},$$
(43)
where $`N_v`$ is the total number of vortices. We then obtain from Eqs. (39) and (40)

$$B_𝐆=N_v\mathrm{\Phi }_0\left[\frac{\text{Det}Q+Q_{\alpha \beta }ϵ^{\alpha \gamma }ϵ^{\beta \delta }G_\gamma G_\delta }{G^4+G^2Q_{\alpha \beta }\delta _{\alpha \beta }+\text{Det}Q}\right],$$
(44)
with $`ϵ^{12}=1=ϵ^{21}`$, $`ϵ^{11}=0=ϵ^{22}`$, and $`\mathrm{\Phi }_0=hc/2e`$ is the quantum of flux. Therefore, the free energy for a vortex lattice per unit volume becomes

$$=\frac{1}{8\pi }(\mathrm{\Phi }_0n_v)^2\underset{𝐆}{}\left[\frac{\text{Det}Q+Q_{\alpha \beta }ϵ^{\alpha \gamma }ϵ^{\beta \delta }G_\gamma G_\delta }{G^4+G^2Q_{\alpha \beta }\delta _{\alpha \beta }+\text{Det}Q}\right]^2\left[1+\frac{G_\alpha Q_{\alpha \beta }G_\beta }{\text{Det}Q}\right].$$
(45)
This has an approximate but much simpler form as

$$\frac{1}{8\pi }(\mathrm{\Phi }_0n_v)^2\underset{𝐆}{}\frac{G_x^2Q_{yy}+G_y^2Q_{xx}G_xG_y(Q_{xy}+Q_{yx})}{G^4},$$
(46)
which is essentially important for determining the ground state structure of the vortex lattice. This form is exact when $`|\frac{\mathrm{}}{2}(\theta )_𝐆||eA_𝐆|/c`$ which is true. The component $`G=0`$ will give free energy $`\overline{}`$ for average magnetic induction $`B`$. For determination of vortex lattice structure, one should in principle minimize Gibbs free energy $`𝒢=BH/4\pi `$. Here $`H`$ is the applied magnetic field. Beyond $`H_{c_1}`$, magnetic field penetrates the superconductor almost fully. Thus $`BH`$, specially so in high $`T_c`$ superconductors, since $`H_{c_1}H_{c_2}`$. $`B`$ does not vary much for different vortex lattice structures for a given $`H`$ as we see in our numerical study that the ratio $`(HB)/B10^7`$. We, therefore, minimize $`\mathrm{\Delta }=\overline{}`$ ie. that part of the free energy which depends on $`G`$, in Eq. (46) for different choices of nonzero $`G`$’s corresponding to different structures and with a cutoff $`G\pi /\xi `$.

## III Numerical Study

The values of the phenomenological parameters that we have used for the numerical computation of $`Q_{\alpha \beta }(q_x,q_y)`$ and later for the free energy are taken from angle resolved photoemission experiments , penetration depth measurement  and band structure calculations  for high-Tc compounds. These are as follows: $`t=1150`$K, $`t^{}=0.48t`$, $`\mathrm{\Delta }_0=400`$K, $`a=3.8\AA `$, $`d=10\AA `$, $`\xi =20\AA `$, $`\lambda =1600\AA `$, and $`\mu =1.33t`$ which corresponds to doping concentration $`x0.19`$.

Using standard Gaussian quadrature, we integrate over $`k_x`$ and $`k_y`$ in Eq. (16) to obtain $`\chi _{\alpha \beta }^\mathrm{p}(𝐪)`$. In Fig. 5 we have shown the dependence of paramagnetic susceptibilities (a) $`\chi _{xx}^\mathrm{p}`$, (b) $`\chi _{yy}^\mathrm{p}`$, and (c) $`\chi _{xy}^\mathrm{p}`$ in units of $`\chi ^d`$ for positive $`q_x`$ at different positive values of $`q_y`$. Susceptibilities for negative values of $`q_x`$ and $`q_y`$ can be obtained by using the symmetries discussed in the previous section. We see that $`\chi _{xx}^\mathrm{p}(q_x,q_y)\chi _{yy}^\mathrm{p}(q_x,q_y)`$ in general. This strong anisotropy in the diagonal susceptibilities is due to the strong $`𝐤`$ dependence of nature of $`\mathrm{\Delta }_𝐤`$ and the $`𝐤`$-dependent $`\gamma _{\alpha \beta }(𝐤)`$ (18)–(21). The diagonal susceptibilities are large compared to off-diagonal one.

We numerically fit, guided by the semianalytical form in Eqs.(28)–(34), to obtain the approximate functional form of $`\chi _{\alpha \beta }^\mathrm{p}(q_x,q_y)`$ for $`q_xa,q_ya0.3`$ as

$`\chi _{xx}^\mathrm{p}(q_x,q_y)`$ $`=`$ $`\{\begin{array}{ccc}\gamma \left[0.31(\delta |q_x|a)+0.14(\delta q_xa)^20.35(\delta q_xa)^2\mathrm{ln}|\delta q_xa|\right]\hfill & \mathrm{if}\hfill & |q_x||q_y|,\hfill \\ \gamma [0.35(\delta |q_y|a)0.14(\delta q_ya)^2+0.10(\delta q_ya)^2\mathrm{ln}|\delta q_ya|\hfill & & \\ +(0.10+\frac{0.21}{\delta |q_y|a})(\delta q_xa)^2+(0.16+\frac{0.07}{\delta |q_y|a})(\delta q_xa)^2\mathrm{ln}|\delta q_xa|]\hfill & \mathrm{if}\hfill & |q_x||q_y|,\hfill \end{array}`$ (51)
$`\chi _{yy}^\mathrm{p}(q_x,q_y)`$ $`=`$ $`\chi _{xx}^\mathrm{p}(q_y,q_x),`$ (52)
$`\chi _{xy}^\mathrm{p}(q_x,q_y)`$ $`=`$ $`\gamma \left[(0.11+{\displaystyle \frac{0.02}{q_>}})q_<+(0.15{\displaystyle \frac{0.18}{q_>}})q_<^2+(0.07{\displaystyle \frac{0.12}{q_>}})q_<^2\mathrm{ln}q_<\right]\mathrm{sign}(q_xq_y)`$ (53)
with $`q_{>,<}=\mathrm{max},\mathrm{min}(|q_x|,|q_y|)\delta a`$, $`\gamma =\frac{\lambda ^2}{d}(\frac{4\pi e^2}{c^2\mathrm{}^2})t`$ whose numerical value is 1.18, and the parameter $`\delta =t/\mathrm{\Delta }_0`$. These phenomenological forms can be explained from the semianalytical expressions (28)–(34) as follows. (i) Firstly, why do $`\chi _{xx}^\mathrm{p}(q_x,q_y)`$ and $`\chi _{yy}^\mathrm{p}(q_x,q_y)`$ not depend on the signs of $`q_x`$ and $`q_y`$, and $`\chi _{xy}^\mathrm{p}(q_x,q_y)`$ does depend on sign$`(q_xq_y)`$? This is due to the symmetry reason discussed following Eq.(22). (ii) Why does $`\chi _{xx}^\mathrm{p}(q_x,q_y)`$ does mainly depend on whether $`|q_x||q_y|`$ or not? This can be understood by the following exercise. We find a term from Eq.(28) as $`|q_y+q_x|`$ assuming $`q_x,q_y0`$. The corresponding term will be $`|q_yq_x|`$ when we consider the contribution from 2nd quadrant. When we add these two contributions, we see that the sum depends on the greater of $`q_x`$ and $`q_y`$. (iii) Following the argument above in (ii), the difference between the two terms is smaller of $`q_x`$ and $`q_y`$. This is the reason why $`\chi _{xy}^\mathrm{p}(q_x,q_y)`$ depends mainly on the smaller of $`|q_x|`$ and $`|q_y|`$. (iv), Since $`\chi _{xx}^{\mathrm{p},1}(q_x,q_y)\chi _{yy}^{\mathrm{p},1}(q_x,q_y)`$, and $`\chi _{xx}^\mathrm{p}(q_x,q_y)=\chi _{yy}^\mathrm{p}(q_y,q_x)`$ for symmetry reasons, the dependence of $`\chi _{xx}^\mathrm{p}(q_x,q_y)`$ on $`q_x`$ and $`q_y`$ is asymmetric. (v) The linear, quadratic, and the logarithmic dependences on $`q`$ follow from in the expressions (28)–(34).

We next numerically perform angular integrals in equations (28)–(38) along with the contributions from other three quadrants to obtain semianalytical susceptibilities and then compare with the fully numerically obtained susceptibilities. In Fig. 6 we have shown $`\chi _{xx}^\mathrm{p}(q_x,0)`$ and $`\chi _{yy}^\mathrm{p}(q_x,0)`$ evaluated in the two ways. The linear approximation of energies in the analytical expressions is a good approximation for determining linear dependence on $`q_x`$ as we see in Fig. 6 that they agree for very low $`q_x`$. They however differ for higher $`q_x`$ since our analytical expressions are not consistent in determining quadratic dependences on $`q`$ as we have neglected higher order $`k`$ dependences to quasiparticle energy. It is however clear that $`\chi _{xx}^\mathrm{p}\chi _{yy}^\mathrm{p}`$ which is our main result.

We consider a face centered rectangular vortex lattice (as shown in Fig. 2) with area of the unit cell $`\stackrel{~}{A}=2\mathrm{\Phi }_0/B`$, in general. Angle $`\theta `$ determines the sides of the rectangle with a fixed area. The sides of the rectangle are $`a_1=[\stackrel{~}{A}\mathrm{tan}\theta ]^{1/2}`$ and $`a_2=[\stackrel{~}{A}/\mathrm{tan}\theta ]^{1/2}`$. We then readily obtain reciprocal lattice vectors for a vortex lattice, in general, to be

$$G_{mn}(B,\theta )=(n+m)\frac{2\pi }{a_1}\widehat{e}_x+(nm)\frac{2\pi }{a_2}\widehat{e}_y,$$
(54)
where $`n`$ and $`m`$ are integers (both positive and negative) including zero. If the vortex lattice makes an angle $`\varphi `$ with the underlying atomic lattice, we find

$`G_{mn}(B,\theta ,\varphi )`$ $`=`$ $`\widehat{e}_x\left[(n+m){\displaystyle \frac{2\pi }{a_1}}\mathrm{cos}\varphi (nm){\displaystyle \frac{2\pi }{a_2}}\mathrm{sin}\varphi \right]`$ (56)
$`+\widehat{e}_y\left[(n+m){\displaystyle \frac{2\pi }{a_1}}\mathrm{sin}\varphi +(nm){\displaystyle \frac{2\pi }{a_2}}\mathrm{cos}\varphi \right].`$
The lattice is a centered square for $`\theta =45^o`$ and triangular when $`\theta =60^o`$. There is symmetry of rotation about $`\varphi =45^o`$, since the lattice is considered as centered rectangular. We therefore need to determine free energy for $`45^o\theta <90^o`$ and $`0\varphi 45^o`$.

We then numerically compute the free energy per vortex (without the single vortex energy which does not depend on structure), $`\mathrm{\Delta }F=d\mathrm{\Delta }/n_v`$ using Eqs.(46) and (56) as functions of the parameters $`\theta `$, $`\varphi `$, and $`B`$. The reciprocal lattice vector $`𝐆`$ changes with the change of any one or more of the parameters. Thus $`\mathrm{\Delta }F`$ differs for different combination of these parameters $`(\varphi ,\theta ,H)`$. In Fig. 7 we have shown the dependence of $`\mathrm{\Delta }F`$ at a low field $`B=2`$ Tesla as a function of $`\varphi `$ for the angles $`\theta =60^o`$ and two neighboring angles $`\theta =58^o`$ and $`62^o`$ (on either side of $`\theta =60^o`$). It is clear that $`\mathrm{\Delta }F`$ is a minimum for the triangular lattice. We notice also that $`\mathrm{\Delta }F`$ is a minimum for the triangular lattice when $`\varphi =0^o`$ and $`30^o`$ which in fact correspond to the same lattice configuration. We likewise find that in the whole of the low field regime, the ground state configuration of the vortex lattice is triangular with one of its arms parallel to one of the crystal axes.

Interestingly, the orientation of the lattice changes discontinuously as we increase the magnetic field though the structure continues to be triangular. In Fig. 8 we have shown the dependence of $`\mathrm{\Delta }F`$ on $`\varphi `$ for the triangular lattice configuration at three chosen fields 2, 5, and 8 Tesla. At nearly about 5 Tesla field, $`\mathrm{\Delta }F`$ is minimum for all $`0^o`$, $`30^o`$, $`15^o`$, and $`45^o`$ orientations; the latter two angles correspond to same lattice configuration, like the former two angles. On the other hand at the field of 8 Tesla, $`\mathrm{\Delta }F`$ is minimum for $`\varphi =15^o`$ and $`45^o`$ only. The triangular vortex lattice changes its orientation discontinuously at about 5 Tesla field. While the triangular lattice has one of its arms parallel to one of the crystal axes at lower field, it aligns to one of the crystal axes by $`45^o`$ at higher field. We understand this discontinuous transition by comparing the energies contributed to $`\mathrm{\Delta }F`$ by the $`𝐆`$ vectors of the lowest magnitude (since they contribute most to the free energy) for these two preferred orientations. Considering the symmetries of the susceptibilities, it is sufficient that we consider only those $`𝐆`$ vectors which have positive $`G_x`$. We thus consider 3 $`𝐆`$ vectors for each of these two orientations. These are (a) $`(1/2,\pm \sqrt{3}/2)G`$ and $`(1,0)G`$ for $`\varphi =0^o`$ and (b) $`\frac{1}{2\sqrt{2}}(\sqrt{3}1,\sqrt{3}+1)G`$, $`\frac{1}{2\sqrt{2}}(\sqrt{3}+1,\sqrt{3}1)G`$, and $`(1/\sqrt{2},1/\sqrt{2})G`$ for $`\varphi =45^o`$, where the length of the smallest $`𝐆`$ vector $`G=2\pi (2/\sqrt{3})^{1/2}(B/\mathrm{\Phi }_0)^{1/2}`$. In Fig. 9, we have shown the energy E contributed by these individual $`𝐆`$ vectors to $`\mathrm{\Delta }F`$ for 2, 5, and 8 Tesla fields. We find the total energy contributed by above 3 $`𝐆`$ vectors for $`\varphi =0^o`$ and $`45^o`$ orientations as (a) 296.05K and 296.51K for $`B=2`$ Tesla, (b) 282.96K and 283.11K for $`B=5`$ Tesla, and (c) 274.15K and 273.60K for $`B=8`$ Tesla respectively. Clearly, the triangular lattice makes an orientational transition at about 5 Tesla field.

To understand the field scale 5 Tesla for the above orientational transition, we compare the energy contributed by the above 3 $`𝐆`$ vectors for each of the preferred orientations. The ratio of these energies can be expressed as a function of $`\alpha =tGa/\mathrm{\Delta }_0`$ using equations (46)–(56). This is given by

$$\frac{E_1}{E_2}=\frac{3f_1(\alpha )}{3f_2(\alpha )}$$
(57)
where $`E_1(E_2)`$ is the energy contributed by the above corresponding 3 $`𝐆`$ vectors of triangular lattice with $`0^o(45^o)`$ orientation. In Fig. 10, we have shown the ratio $`f_1(\alpha )/f_2(\alpha )`$ as a function of $`\alpha `$. The orientational transition takes place when the ratio is unity. This corresponds to $`\alpha _c=tG_ca/\mathrm{\Delta }_00.37`$. Therefore the critical field at which the transition take place, $`B_1(0.37/2\sqrt{2}\pi )^2\sqrt{3}(\mathrm{\Delta }_0/at)^2\mathrm{\Phi }_05.2T`$.

The structure of the vortex lattice remains triangular with $`45^o`$ orientation to the crystal lattice as we have shown $`\mathrm{\Delta }F`$ in Fig. 11 for a field as high as $`25`$ Tesla. However, it makes yet another discontinuous structural transition to a centered square lattice with its axes parallel to the crystal axes at yet another critical field $`B_2`$ whose value is about 28 tesla. Figure 12 shows that $`\mathrm{\Delta }F`$ is minimum for $`\theta =45^o`$ and $`\varphi =0^o`$ at $`B=28`$ Tesla. The overall phase diagram for the ground state of the vortex lattice structure at $`T=0`$ is shown in figure 1.

## IV Concluding Remarks

We conclude by briefly discussing a number of questions such as the nature of the approximations used, the effect of nonzero temperature, consequences of the transitions, their observability, the reason why our conclusions differ from those found earlier, and the structural transitions experimentally observed.

We have calculated the ground state energy assuming effectively that the interaction between two vortices is unaffected by the presence of other vortices. This is obviously a low vortex density approximation which seems quite reasonable since the dimensionless ratio $`(n_v/n)`$ is about 1/2500 for a field of a Tesla. However, we have not calculated the higher order corrections which while nominally of higher order in $`(n_v/n)`$ might have large or even divergent coefficients. Since the vortex vortex interaction depends on quasiparticle-quasihole susceptibility, a change in their spectrum due to the supercurrent (the Volovik effect  ) could have serious consequences. Here, we would like to make two points. Firstly, a calculation by Amin, Affleck and Franz , using a semiclassical approximation to include the nonlinear effect of the magnetic field a la Volovik, finds that this has little effect on the structural transformations calculated by them. Secondly, all the recent fully quantum mechanical calculations  of the density of Dirac quasiparticle states in a vortex lattice in the London limit find that for Bravais lattices, the density of states vanishes linearly with energy as in the absence of a magnetic field; only the slope changes. A general argument for this has been presented by Vishwanath . For these two reasons, we believe that our low density approximation is reliable.

In the London approximation, the vortex cores are treated as $`\delta `$ functions. In reality, they have a width of order the coherence length. We believe that the consequences of this approximation, at least for the low field structural transition, are small. The reason is that the structure sensitive part of the energy arises from the difference in contribution of the smallest reciprocal lattice vectors (Section III). For these, $`G\xi (1/10)`$ at typical magnetic fields, so that the phenomenological assumption of a Gaussian vortex core with width $`\xi `$ will make a negligible difference to the structure sensitive part.

We have calculated only the ground state energy of the vortex lattice in this paper. At any nonzero temperature, there are obviously entropic contributions which could change the magnetic field at which the structural transition occurs, as a function of temperature. Here, we note that since both the structures (below and above 5T) are identical (triangular), and the structure difference sensitive part of the energy is a tiny fraction $`(<(1/1000)^2)`$ of the vortex interaction energy, the elastic fluctuations in both structures are expected to be identical to order $`(1/1000)^2`$ so that the transition field should not be affected by temperature, so long as the input parameters (e.g., $`\xi _𝐤,\mathrm{\Delta }_𝐤,\lambda `$ ) do not change with $`T`$. The same cannot be said of the high field $`(28T)`$ triangular to centred square lattice transition, because one has a tight packed structure and the other not. The expectation is that the former has lesser elastic fluctuations than the latter, so that the transition field boundary should shift to lower values with increasing temperature. However, this conclusion is tempered by the fact that the London approximation is unreliable at these high fields when vortex cores get close to each other so that our basic result may not be that reliable.

One interesting consequence of the orientational transition at 5T, which might be measurable, is the change in the very low energy density of quasiparticle states. At least for a square lattice, Vishwanath  has shown that there are quasiparticle states with linear dispersion and that there is a very small gap arising from higher order terms in the quasiparticle velocity. If this kind of result carries through for a triangular lattice, then it might be an experimental way of observing the transition.

We have calculated here, for the first time, the actual energy of the structure sensitive term (Section III) and have found it to be small, of order a few degrees per vortex. Because of this reason, the transition might be difficult to observe, since pinning energies of larger size are generally present , unless the system is extremely perfect.

We have discussed in detail (Section II and III) the reason why our results differ from those obtained earlier. Basically, it has to do with the anisotropy of the nonlocal current current susceptibility, i.e. the fact that $`\chi _{xx}^\mathrm{p}(𝐆)\chi _{yy}^\mathrm{p}(𝐆)`$. The reason for this essentially is that we have an anisotropic superconductor. The $`\chi _{xx}^\mathrm{p}`$ and $`\chi _{yy}^\mathrm{p}`$ functions are plotted in Fig.5.

The question of a non triangular structure of the vortex lattice in cuprates has attracted considerable experimental attention , especially since it has become established that they are $`d_{x^2y^2}`$ superconductors. Earlier small angle neutron scattering measurements  were on highly twinned 123 crystals, so that the observation of four fold diffraction symmetry does not imply rectanglar /square lattice. Moreover, the positional order is very poor. A more recent experiment  on untwinned 123 single crystals shows much better translational order (higher $`G`$ peaks are resolved for the first time) and a triangular lattice with axes oriented along a, distorted because of $`a`$-$`b`$ asymmetry. The authors find no structural transitions upto 4T with field along $`c`$-axis. They however find a transition from a triangular lattice oriented along $`a`$ to one oriented along the $`b`$ axis at a field of about 3.8T, oriented at 33o to the $`c`$ axis. This is certainly quite different from the transition to a triangular lattice at 45o to $`a`$ axis at 5T predicted by us. As Johnson et al point out  the observed transition could be due to the presence of chains in 123, and the novel $`ab`$-$`c`$ anisotropy caused by it (which may have a strong effect on many physical properties). In order to seriously explore our conclusions, one needs to do experiments on cuprates without chains and ideally with tetragonal symmetry, as again mentioned by Johnson et al .

## ACKNOWLEDGMENTS

We acknowledge support from the JNC India and the US-India ONR funded project No. N00014-97-0988.