format: text

text:
# Numerical Study of Wave Propagation in Uniaxially Anisotropic Lorentzian Backward Wave Slabs

## I Introduction

Metamaterials have received much attention during the last years, because they possess unusual electromagnetic properties, like, for example, the opposite directions of phase and group velocities. Double negative (DNG) materials have negative permittivity and permeability, and they belong to the class of metamaterials. These media that are capable of supporting backward waves, have been also called backward wave (BW) media in the literature Lindell et al. (2001). In BW media, the refraction phenomenon is anomalous in the sense that the power flow is refracted negatively, i.e. to the same side of the normal of the interface. As discussed in Lindell et al. (2001), it is not necessary for the medium to be a DNG medium to be able to support backward waves, because anomalous refraction can also be realized with anisotropic media with only one negative material parameter. We will use the term BW medium throughout the rest of the paper.

The pioneering work on BW materials by Veselago Veselago (1968), where slab lenses were mentioned, has gained much attention during recent years, despite some differing opinions on the subject Valanju and Walser (2002). Isotropic BW materials are often called Veselago materials. The possibility to realize a perfect lens with isotropic BW slabs was discussed by Pendry in Pendry (2000). Numerical and theoretical considerations of wave propagation in isotropic BW slabs excited with a line current above the slab were presented by Ziolkowski and Heyman in Ziolkowski and Heyman (2001). Guidance of waves in a slab of uniaxially anisotropic metamaterial has been theoretically discussed by Lindell and Ilvonen in Lindell and Ilvonen (2002).

Controversial opinions regarding negative refraction effect and the perfect lens call for more detailed studies of the wave phenomena in backward-wave media. In this paper, wave propagation through uniaxially anisotropic BW slabs is numerically studied, and comparison is made with the theory presented in Lindell et al. (2001). The theory predicts that there are regions in certain BW media, where the wave vector becomes complex, thus resulting in exponentially decaying waves. These regions are bounded by the asymptotes of the wave vector surfaces, which can be shown to be hyperbola. We study these phenomena numerically using the finite-difference time-domain (FDTD) method in a 2D-problem of a line current radiating in the vicinity of a BW slab. Also, the existence of surface waves on the interface between free space and BW medium is demonstrated with an example case. The BW medium is realized with Lorentzian constitutive parameters having a single pole pair.

An example problem with some theoretical discussion are presented in section $`3`$. Results from the numerical simulations are shown and discussed in section $`4`$. Our numerical simulations show that the wave propagation and refraction phenomena heavily depend on the parameter choices of the BW medium and are qualitatively in agreement with the theory.

## II Numerical Model of Dispersive Medium

The constitutive relations for a frequency dispersive isotropic medium read

$$𝐃=ϵ(\omega )𝐄,𝐁=\mu (\omega )𝐇.$$
(1)
Negative permittivity and permeability are realized with the Lorentz medium model. The expressions for the permittivity and permeability are of the form

$`ϵ(\omega )`$ $`=`$ $`ϵ_0\left(1+{\displaystyle \frac{\omega _{pe}^2}{\omega _{0e}^2\omega ^2+j\mathrm{\Gamma }_e\omega }}\right),`$
$`\mu (\omega )`$ $`=`$ $`\mu _0\left(1+{\displaystyle \frac{\omega _{pm}^2}{\omega _{0m}^2\omega ^2+j\mathrm{\Gamma }_m\omega }}\right).`$ (2)
This model corresponds to a realization of BW materials as mixtures of conductive spirals or omega particles, as discussed in Tretyakov et al. (2002). In this artificial material both electric and magnetic polarizations are due to currents induced on particles of only one shape, which provides a possibility to realize the same dispersion rule for both material parameters, as in (2). Note that the medium realized by Smith et. al. is built using different principles Smith et al. (2000). For the uniaxial materials that we consider in this paper we assume that the negative components of the material parameters are realized by small uniaxial spiral inclusions (racemic arrays with equal number of right- and left-handed particles) and possess frequency dispersion defined by (2). The positive components of the material parameters are equal to the free-space permittivity and permeability values, assuming that there are no particles oriented along these axes.

Equations (1) and (2) form the basis of the used FDTD model for BW materials. The two most important known FDTD methods for modeling dispersive materials like Lorentz materials are the recursive convolution method and the auxiliary differential equation method, where the constitutive parameters are expressed with the help of susceptibility. In the first method, $`𝐃`$ and $`𝐄`$ and $`𝐁`$ and $`𝐇`$ are related through a convolution integral. This approach is rather tedious. Another possibility is to use the auxiliary differential equation technique, which is slightly easier to implement. In this latter method, the polarization current associated to each Lorentz pole pair is introduced. These two models are discussed in detail in Taflove and Hagness (2000).

A third method to discretize fields in Lorentz medium, classified as direct integration method in Young and Nelson (2001), is based on the direct discretization of the PDE representing the time-domain equivalent of the simplified frequency-domain constitutive relation. The proposed discretization scheme is a modification of this method. The idea is to transform (1) into the time domain using the relation $`j\omega /t`$ with one integration before discretization using center differences. Usually, FDTD models based on the constitutive relation are directly (after multiplication with the denominator) discretized, as discussed in a summary of FDTD algorithms for dispersive media in Young and Nelson (2001). We found in K$`\ddot{a}`$rkk$`\ddot{a}`$inen and Maslovski (2003) that one integration prior to discretization leads to much better accuracy and, what is also important, to considerably better stability properties than the direct discretization without integration. In the numerical simulations, we have used the model which is discussed in detail in K$`\ddot{a}`$rkk$`\ddot{a}`$inen and Maslovski (2003).

## III An Example Problem and Theoretical Discussion

Consider a $`z`$-directed line current in free space located at a distance $`d_s`$ from a BW slab of thickness $`d`$. Let the interface between free space and the BW slab be located at $`y=0`$. The problem space is two-dimensional, with the field components $`H_x,H_y`$, and $`E_z`$. The peak of the incident spectrum is at $`\omega _p=5.010^9`$ rad/s, and the parameters in (2) are the following: $`\omega _{0e}=\omega _{0m}=1.010^9`$ rad/s, $`\omega _{pe}^2=\omega _{pm}^2=4.810^{19}`$ (rad/s$`)^2`$, $`\mathrm{\Gamma }_e=\mathrm{\Gamma }_m=0`$. With these choices, we obtain $`ϵ(\omega )=\mu (\omega )`$ for all $`\omega `$ and $`ϵ(\omega )/ϵ_0=\mu (\omega )/\mu _0=1`$ at $`\omega =\omega _p`$. The spatial resolution $`\mathrm{\Delta }x=\mathrm{\Delta }y=1.5`$ cm is used throughout the simulations. To be able to demonstrate the properties of BW materials, the incident spectrum is quite narrow, so that the relative constitutive parameters are close to minus one for the frequencies having significant spectral content. Absorbing boundary conditions are used to terminate the computational domain at the outer boundaries of the lattice. For simplicity, we have used Liao’s third order ABC, although more sophisticated ABC’s are available. The use of usual ABC’s requires a small gap between the outer boundary of the computational space and the BW material slab. The chosen coordinate axes and the problem geometry is shown in Figure 1.

Next, we briefly present some theoretical results from Lindell et al. (2001) that are important here for comparison purposes with our numerical results. Wave propagation in BW slabs is studied with different value combinations of the medium parameters $`\mu _x`$, $`\mu _y`$ and $`ϵ_z`$. For our TE-polarized case, the Poynting vector $`𝐒_{TE}`$ can be shown to read Lindell et al. (2001)

$$𝐒_{TE}=\frac{\overline{\overline{\mu }}𝐤_{TE}}{2k_0\eta _0\mu _0\mu _x\mu _y}|E_0|^2,$$
(3)
where $`|E_0|`$ is the amplitude of the TE-polarized electric field of the plane wave, $`𝐤_{TE}`$ is the wave vector with two cartesian components $`k_x`$ and $`k_y`$, $`k_0=\omega \sqrt{ϵ_0\mu _0}`$ is the free space wave number, and $`\eta _0=\sqrt{\mu _0/ϵ_0}`$ is the free space wave impedance. Denoting the angle between the outwards-pointing unit normal vector $`𝐮_y`$ and the wave vector $`𝐤_{TE}`$ by $`\theta `$, one can derive the dispersion equation and solve it for the wave number as in Lindell et al. (2001). The result is

$$k_{TE}(\theta )=k_0\sqrt{\frac{\mu _xϵ_z}{\mathrm{cos}^2\theta +\frac{\mu _x}{\mu _y}\mathrm{sin}^2\theta }}.$$
(4)
For any given set of parameters, we may plot the projection curves of the wave vector surfaces in $`xy`$-plane. For certain choices of the medium parameters, the wave number becomes complex, resulting in exponentially decaying waves inside the BW material.

Any physical power flow must be directed downwards away from the source. This requires that

$$𝐮_y𝐒_{TE}<0.$$
(5)
The necessary condition for any transmission is

$$\frac{𝐮_y𝐤_{TE}}{\mu _x}<0\frac{k_y}{\mu _x}<0.$$
(6)
Clearly, the phase velocity is directed oppositely to the power flow provided that $`\mu _x<0`$. Negative refraction of the Poynting vector requires that

$$k_x𝐮_x𝐒_{TE}<0\mu _y<0.$$
(7)
For the interface problem, we will also need to know conditions for the existence of surface waves. The input impedance on the surface filled by a uniaxial material is, for TE-polarized fields,

$$Z_{\mathrm{inp}}=\frac{\omega \mu _x}{\beta _{TE}},$$
(8)
where the normal component of the propagation factor reads

$$\beta _{TE}=\sqrt{\frac{\mu _x}{\mu _y}\left(\omega ^2ϵ_z\mu _yk_x^2\right)}$$
(9)
with the square root branch defined by $`\mathrm{Im}(\beta _{TE})<0`$. $`k_x`$ is the propagation factor along the surface. Thus, the eigensolutions for an interface between this medium and free space satisfy the following equation:

$$\frac{\omega \mu _x}{\beta _{TE}}+\frac{\omega \mu _0}{\beta _0}=0$$
(10)
where

$$\beta _0=\sqrt{\omega ^2ϵ_0\mu _0k_x^2},$$
(11)
with $`\mathrm{Im}(\beta _0)<0`$.

Let us consider now conditions for the existence of surface waves along this interface. In this case both betas must be imaginary, of course with negative imaginary parts: $`\beta _{TE}=j\alpha _{TE}`$ and $`\beta _0=j\alpha _0`$, where $`\alpha _{TE}>0`$ and $`\alpha _0>0`$. The eigenvalue equation becomes

$$\frac{\omega \mu _x}{\alpha _{TE}}+\frac{\omega \mu _0}{\alpha _0}=0.$$
(12)
Obviously, if all the material parameters are positive, this equation has no solutions, but if $`\mu _x<0`$, surface wave solutions are possible. This is well known for interfaces with free-electron plasma.

Five different cases will be considered in the following. In the first case, we choose $`\mu _x>0`$, $`\mu _y<0`$, and $`ϵ_z<0`$. The wave vector surface is a two-sheeted hyperboloid with axis parallel to the $`x`$-axis. The asymptotes of the hyperbola in $`xy`$-plane divide the plane into regions of complex and real wave vectors. Waves that are propagating parallel to $`y`$-axis are supposed to decay exponentially inside the slab.

In the second case, $`\mu _x>0`$, $`\mu _y<0`$, $`ϵ_z>0`$. The wave vector surface is a two-sheeted hyperboloid with the axis parallel to the surface normal. In our 2D-case, real wave vectors exist inside the region bounded by the asymptotes of the hyperbola associated to the wave vector surfaces.

As a third case, we consider the situation complementary to the second case in the sense that the signs of all the parameters are changed. Notice that this does not affect the shape of the wave vector curves. We have $`\mu _x<0`$, $`\mu _y>0`$, and $`ϵ_z<0`$. In fact, the Poynting vector is refracted positively in this case (see (7)), but this is an interesting case anyway because of the aforementioned contrast with respect to the second case.

The fourth case consists of the usual isotropic BW slab ($`\mu _x=\mu _y<0`$, and $`ϵ_z<0`$) where focusing and negative refraction phenomena are present. The wave vector curves are ellipses (in our case of two equal parameters they are circles), and real wave vectors exist everywhere inside the slab.

The fifth case is specially chosen to show the existence of surface waves in the case when $`\mu _x<0`$, $`\mu _y<0`$, and $`ϵ_z>0`$. We can easily see from (4) that there are neither real wave vectors nor backward waves. However, we have found that surface waves on the interface are easily excited in this case. Let us now present the numerical results.

## IV Numerical Results and Comparison with the Theory

In the first four cases, we show the electric field distribution at three suitable chosen increasing time steps to illustrate the wave propagation and refraction phenomena. In the fifth case, we illuminate a rectangular cylinder to see the surface waves. Whenever a constitutive parameter is said to be positive, it is supposed to be a constant and equal to the free space permittivity or permeability. Negative material parameters obey the Lorentzian dispersion rule and equal $`ϵ_0`$ or $`\mu _0`$ at the center frequency.

### IV.1 Case I: $`\mu _x>0`$, $`\mu _y<0`$, $`ϵ_z<0`$

In this case, the theory shows that the wave vectors are complex inside the slab within a region bounded by the asymptotes of a hyperbola. Inspection of Figure 2 c) reveals that there is indeed a region in the slab, where the electric field in negligible all the time. There are some fields within the slab just under the source. The hyperbola-shaped wavefronts propagate obliquely downwards inside the slab and the power flow is refracted negatively. The distance of the source from the first interface is discretized with $`10`$ cells, and the thickness of the slab corresponds to $`80`$ cells.

### IV.2 Case II: $`\mu _x>0`$, $`\mu _y<0`$, $`ϵ_z>0`$

In Figure 3 a), a cylindrical wave is penetrating into the BW slab. In Figure 3 b), some numerical dispersion is visible. There are significant fields inside the region where the theory yields real wave vectors, while the fields are rather small elsewhere. Despite some dispersive effects, the wavefronts of constant field value are reminiscent of hyperbolas. The phase velocity is directed downwards. Some weak focusing of the power flow is seen in Figure 3 b). The phase velocity inside the slab is directed downwards, as can be seen from the theory.

### IV.3 Case III: $`\mu _x<0`$, $`\mu _y>0`$, $`ϵ_z<0`$

To complete the analysis, we change the signs of the parameters of the second case. Clearly, the wave vector surfaces as defined by (4) are not changed. In fact, the Poynting vector is refracted positively in this case. However, the phase velocity is directed upwards \[see (6) and (7)\]. This phenomenon is clearly seen during the simulation. From Figure 4 we see that the electric field distributions inside the slab are quite similar to those of Figure 3 except that the wavefronts are less distorted in Figure 4. In Figure 4 c), the wavefronts behind the slab are seen to be oblate ellipses with the center on the lower interface of the slab.

### IV.4 Case IV: an Isotropic Slab with $`\mu _x<0`$, $`\mu _y<0`$, $`ϵ_z<0`$

Here we consider the usual isotropic BW (or double negative) slab with all the relative material parameters close to minus one. This case has been studied, for example, by Ziolkowski and Heyman in Ziolkowski and Heyman (2001) in the case of Drude slabs. We obtained quite similar results with our alternative discretization technique in the case of isotropic Lorentz medium. The electric field distributions are shown in Figure 5. We can calculate the positions of the foci from the slab thickness $`d`$ and the distance $`d_s`$ of the source from the interface. Notice that we must have $`d_s<d`$ to have a focus inside the slab. The foci should appear at $`y=d_s`$ inside the slab and at $`y=(2dd_s)`$ behind the slab. The appropriate derivations can be found in Ziolkowski and Heyman (2001).

From Figure 5 a) we see that the electric field is concentrated in the expected position of the first focus. In Figure 5 c), the second focus behind the slab is also visible. It takes some time for the foci to develop, as is seen from Figure 5 b), where the second focus is not yet seen. The incident spectrum has some small components for which the relative material parameters are not exactly minus unity. Hence the wavefronts are not perfect circles as predicted by the theory. Anyway, these results are in agreement with the theory concerning negative refraction of the Poynting vector. However, no “perfect” focusing has been observed, meaning that the focus area is always not smaller that about half wavelength. Steady state solutions for the foci were not obtained. This same observation was also made by Ziolkowski and Heyman in Ziolkowski and Heyman (2001).

### IV.5 Case V: $`\mu _x<0`$, $`\mu _y<0`$, $`ϵ_z>0`$

For this set of material parameters, the theory predicts that the waves decay exponentially everywhere inside the slab. However, we have found that surface waves are easily excited in this case, in accordance with the theoretical prediction, see s̊urface. To see this phenomenon, we illuminate a rectangular cylinder with a pulse having a slightly broader spectrum. The electric field distribution induced on the surface of the cylinder is shown in Figure 6.

Figure 6 a) shows how the surface waves begin to develop. The source is still on. In Figure 6 b), the source above the cylinder has been switched off. The surface waves have propagated to the opposite side of the cylinder as well. In agreement with the theory, there are no fields inside the cylinder except for the immediate vicinity of the surface of the cylinder.

## V Conclusions

Wave propagation and refraction phenomena in uniaxially anisotropic BW slabs have been studied with the FDTD method. Special attention was paid to the shape of the wavefronts and to the regions inside the slabs, where the wave vector becomes complex, thus resulting in exponentially decaying waves. The numerical results for anisotropic BW slabs were seen to qualitatively agree well with the theoretical results. The effects of negative refraction, (imperfect) focusing, and surface wave excitation have been demonstrated. Potentially useful transformations of wave fronts between spherical, elliptical, and hyperbolical can be realized in homogeneous uniaxial backward wave slabs.

## Acknowledgments

Financial support was received from Nokia Foundation and Emil Aaltonen Foundation. Useful discussions with S.A. Tretyakov, S.I. Maslovski and P.A. Belov are gratefully acknowledged.


text:
# Vortex Lattice Structural Transitions: a Ginzburg-Landau Model Approach

## Abstract

We analyze the rhombic to square vortex lattice phase transition in anisotropic superconductors using a variant of Ginzburg-Landau (GL) theory. The mean-field phase diagram is determined to second order in the anisotropy parameter, and shows a reorientation transition of the square vortex lattice with respect to the crystal lattice. We then derive the long-wavelength elastic moduli of the lattices, and use them to show that thermal fluctuations produce a reentrant rhombic to square lattice transition line, similar to recent studies which used a nonlocal London model.

Vortex lattices in superconductors with four-fold anisotropy can exhibit a rhombic to square structural phase transition due to a competition between the anisotropy (favoring a square lattice) and the repulsive vortex-vortex interaction (favoring a triangular lattice). This transition has been observed in the borocarbide family of superconductors $`R`$Ni2B2C (with $`R`$ a rare earth element), whose anisotropy is due to lattice effects eskildsen ; eskildsen2 ; exp3 ; sakata ; vinnikov2 , and more recently in the overdoped cuprate superconductor La1.83Sr0.17CuO4+δ gilardi , whose anisotropy is due to unconventional superconducting order. Recent small angle neutron scattering (SANS) data for LuNi2B2C eskildsen  show the transition line in the $`HT`$ plane curving upward and avoiding $`H_{c2}`$, contrasting with mean-field predictions and underscoring the importance of thermal fluctuations even in this low temperature superconductor. Some features of the data can be understood using a nonlocal London model developed by Kogan and collaborators kogan1 ; kogan2 ; however, the London model is ostensibly only applicable at low magnetic inductions, leaving as an open question the behavior close to $`H_{c2}`$. An alternate explanation for the reentrant behavior of the vortex lattice transition is also available in the literature nakai .

In this Letter we study the effect of thermal fluctuations on the rhombic to square vortex lattice transition using an anisotropic GL model. We first determine the mean-field phase diagram to second order in the anisotropy parameter, and find that even at this level GL theory predicts a novel reorientation transition of the square vortex lattice with respect to the crystal lattice, an effect overlooked in previous GL studies. We then derive the (nondispersive) elastic moduli of the lattices using methods first pioneered by Brandt brandt1  for isotropic superconductors. Finally, these moduli are used to determine the fluctuation renormalized phase boundary in the $`HT`$ plane. For parameters typical of the borocarbides we find a reentrant phase boundary, in agreement with the nonlocal London model results of Ref. gurevich . As GL theory is expected to be valid near $`H_{c2}`$, our results nicely complement and extend the nonlocal London model results.

Accounting for anisotropy, the linearized first GL equation (in dimensionless units) is

$$\mathrm{\Pi }^2\psi \psi ϵ(\mathrm{\Pi }_x^2\mathrm{\Pi }_y^2)^2\psi =0,$$
(1)
where $`𝚷=(\mathbf{}/i\kappa 𝐀)`$ is the covariant derivative and $`ϵ`$ is the dimensionless anisotropy parameter. For the cuprates the anisotropic term can be derived rigorously from the two field $`s`$ and $`d`$-wave GL free energy franz ; heeb  by integrating out the $`s`$-wave term affleck . For the conventional, but highly anisotropic, borocarbides it is a phenomenological term that will be sufficient to capture all the interesting physics. The symmetry breaking properties of the anisotropic term dictate the introduction of the orientation angle $`\varphi `$ of the vortex lattice relative to the underlying crystal lattice—this term provides an effective coupling of the vortex and crystal lattices, and makes structural transitions possible.

We have generalized the virial theorem of superconductivity doria ; klein  to include the anisotropic term in the extended GL theory future , with the result

$`{\displaystyle \frac{1}{4\pi }}𝐇𝐁`$ $`=`$ $`|\psi |^2|\psi |^4+{\displaystyle \frac{𝐁^2(𝐱)}{4\pi }}`$ (2)
$`2ϵRe[\psi ^{}(\mathrm{\Pi }_x^2\mathrm{\Pi }_y^2)^2\psi ],`$
where the brackets $``$ denote spatial averaging.

Using the generalizations of the Abrikosov identities derived with the help of Eq. (2), the free energy can be cast into the compact form

$$F=B^2\frac{(\kappa B)^2}{(2\kappa ^21)\beta _A+1}\left\{1\frac{\mathrm{\Gamma }}{(2\kappa ^21)\beta _A+1}\right\},$$
(3)
where $`\beta _A`$ is the usual Abrikosov parameter and $`\mathrm{\Gamma }`$ is the correction to the isotropic result,

$`\begin{array}{cc}\hfill \mathrm{\Gamma }& =4ϵ\{{\displaystyle \frac{Re|\psi |^2\psi ^{}(\mathrm{\Pi }_x^2\mathrm{\Pi }_y^2)^2\psi }{|\psi |^2^2}}\hfill \\ & +\beta _A(2\kappa ^21){\displaystyle \frac{Re\psi ^{}(\mathrm{\Pi }_x^2\mathrm{\Pi }_y^2)^2\psi }{|\psi |^2}}\}.\hfill \end{array}`$ (4)
Throughout this paper we will use the convenient parameterization $`\rho +i\sigma =(\overline{b}/\overline{a})\mathrm{exp}(i\theta )`$, where $`\overline{b}`$ and $`\overline{a}`$ are the magnitudes of the two basis vectors of the vortex lattice and $`\theta `$ is the apex angle.

Equation (1) was solved perturbatively to second order in the anisotropy parameter $`ϵ`$, extending previous first order results chang1 ; chang2 ; rosenstein . The second order calculation is necessary to justify the use of perturbation theory and at the same time to quantify the accuracy of the first order correction. More importantly, we uncover new and unexpected physics at second order: a reorientation transition of the square vortex lattice with respect to the crystal lattice. The quantities $`\beta _A`$ and $`\mathrm{\Gamma }`$ are given in terms of rapidly converging sums involving $`\rho `$, $`\sigma `$, and $`\varphi `$,

$`\begin{array}{cc}\hfill \beta _A& =\beta _0+(ϵb)\beta _1+(ϵb)^2\left[\gamma _1+\gamma _2+5\beta _1{\displaystyle \frac{3}{4}}\beta _0\right],\hfill \end{array}`$ (5)
$`\begin{array}{cc}\hfill \mathrm{\Gamma }& =8ϵb^2\{\beta _A(2\kappa ^21)[1+3(ϵb)]+\beta _0[1+{\displaystyle \frac{3}{2}}(ϵb)]\hfill \\ & +\beta _1[1+6(ϵb)]+2(ϵb)[\gamma _1+\gamma _2]\},\hfill \end{array}`$ (6)
where the functions $`\beta _0`$, $`\beta _1`$, $`\gamma _1`$, $`\gamma _2`$ are periodic in $`\rho `$ with period $`1`$, and are given by

$`\beta _0`$ $`={\displaystyle \underset{nm}{}^{}}A_{nm},`$ (7)
$`\beta _1`$ $`=Re\left\{e^{4i\varphi }{\displaystyle \underset{nm}{}^{}}A_{nm}\left[8\pi ^2\sigma ^2n^46\pi \sigma n^2+{\displaystyle \frac{3}{8}}\right]\right\},`$ (8)
$`\begin{array}{cc}\hfill \gamma _1& =Re\{e^{8i\varphi }{\displaystyle \underset{nm}{}^{}}A_{nm}[16\pi ^4\sigma ^4n^812\pi ^3\sigma ^3n^6\hfill \\ & +{\displaystyle \frac{3}{4}}\pi ^2\sigma ^2n^4{\displaystyle \frac{45}{16}}\pi \sigma n^2+{\displaystyle \frac{105}{256}}]\},\hfill \end{array}`$ (9)
$`\begin{array}{cc}\hfill \gamma _2& ={\displaystyle \underset{nm}{}^{}}A_{nm}[16\pi ^4\sigma ^4n^4m^412\pi ^3\sigma ^3n^2m^2(n^2+m^2)\hfill \end{array}`$
$`\begin{array}{cc}& +{\displaystyle \frac{3}{4}}\pi ^2\sigma ^2(n^4+m^4+36n^2m^2)\hfill \\ & {\displaystyle \frac{45}{16}}\pi \sigma (n^2+m^2)+{\displaystyle \frac{105}{256}}],\hfill \end{array}`$ (10)
where $`A_{nm}=\sqrt{\sigma }e^{2i\pi \rho (n^2m^2)}e^{2\pi \sigma (n^2+m^2)}`$. The prime denotes a double sum: one over $`n`$, $`m`$ and a second with $`n`$ and $`m`$ replaced by $`(n+1/2)`$ and $`(m+1/2)`$ respectively. We perform a numerical minimization of the free energy (3) with respect to $`\rho `$, $`\sigma `$, and $`\varphi `$ to obtain the mean-field vortex lattice structure as a function of the field and temperature. We find a continuous transition from a rhombic to square phase as $`ϵ`$ is increased. The second order correction only shifts the transition point by about 1% relative to the first order calculation, contrary to previous claims park  which neglected the lifting of the orientation degeneracy by the anisotropic term. Our numerical investigations show unambiguously that a consistent exploration of the relevant phase space and an accurate determination of the global minima of the GL free energy (3) demands the inclusion of the the orientational degree of freedom $`\varphi `$. After the structural phase transition occurs the square vortex lattice is oriented along the  direction, and remains in this preferred orientation until a field $`b_r`$ is reached. At that point another continuous phase transition takes place, with the vortex lattice changing its orientation while retaining its square shape. At a given $`b>b_r`$, two equivalent configurations exist: one with orientation angle $`\varphi >\pi /4`$, and one with $`\varphi ^{}=\pi /2\varphi `$. This reorientation transition occurs only at second order in the anisotropy $`ϵ`$.

Having determined the mean-field phase diagram of the system, we now turn to the elastic properties of the lattices. Assuming that the displacement field $`𝐮`$ is small compared to typical lattice spacings, we calculate the elastic moduli from the GL free energy (3) in the incompressible limit, $`\mathbf{}𝐮=0`$. For simplicity we have focused on straight and parallel vortices; we do not expect the structural transition to significantly affect the tilt modulus $`c_{44}`$. The relevant moduli are essentially non-dispersive brandt1 , allowing us to consider uniform deformations. Now there are many ways to choose the elastic moduli from the independent components of the elastic tensor, and we will follow the lead of Miranović and Kogan mir  and use as the four independent elastic moduli for a rhombic lattice the quantities $`c_{sq}`$, $`c_{66x}`$, $`c_{66y}`$, and $`c_r`$, corresponding to four particular deformations. The so-called “squash” deformation is responsible for the transformation of the rhombic to a square vortex lattice, and the corresponding squash modulus $`c_{sq}`$ vanishes at the transition point. The other relevant deformations are the two pure shears in $`x`$ and $`y`$ directions (which determine $`c_{66x}`$, $`c_{66y}`$) and the rotation of the lattice about the direction of the magnetic field (which determines $`c_r`$).

All elastic moduli are obtained from the second derivative of the GL free energy (3). The elastic energy can be written as

$$E=\frac{1}{2}\mu ^2c(\mu )|_{\mu 0}=\frac{1}{2}\mu ^2_\mu ^2F|_{\mu 0},$$
(11)
where $`c(\mu )`$ is the elastic modulus for each particular deformation. Each elastic modulus is obtained from

$`\begin{array}{cc}\hfill c& =C_{\rho \rho }\rho _\mu ^2+C_{\sigma \sigma }\sigma _\mu ^2+C_{\varphi \varphi }\varphi _\mu ^2\hfill \\ & +2C_{\rho \sigma }\rho _\mu \sigma _\mu +2C_{\rho \varphi }\rho _\mu \varphi _\mu +2C_{\sigma \varphi }\sigma _\mu \varphi _\mu ,\hfill \end{array}`$ (12)
where $`C_{ij}`$ are the terms originating from the derivatives of the GL free energy (3)

$$\begin{array}{c}C_{ij}=(2\delta _{ij})\{\mathrm{\Phi }(b^2F)\beta _{A}^{}{}_{,ij}{}^{}+\mathrm{\Theta }[\mathrm{\Gamma }_{,ij}\mathrm{\Phi }(\beta _{A}^{}{}_{,i}{}^{}\mathrm{\Gamma }_{,j}\hfill \\ \hfill +\beta _{A}^{}{}_{,j}{}^{}\mathrm{\Gamma }_{,i}+\mathrm{\Gamma }\beta _{A}^{}{}_{,ij}{}^{}2\mathrm{\Phi }\mathrm{\Gamma }\beta _{A}^{}{}_{,i}{}^{}\beta _{A}^{}{}_{,j}{}^{})]\}.\end{array}$$
(13)
The derivatives in the right hand side of Eq. (13) are with respect to $`i,j=\rho `$, $`\sigma `$, $`\varphi `$, $`\mathrm{\Phi }=(2\kappa ^21)/[(2\kappa ^21)\beta _A+1]`$, and $`\mathrm{\Theta }=(1b)^2/[(2\kappa ^21)\beta _A+1]^2`$.

The calculated elastic moduli for $`ϵ=0.11`$ and $`\kappa =10`$ (appropriate for the borocarbides) are shown in Fig. 1. The transition from a rhombic to a square vortex lattice is signaled by the vanishing of the squash modulus. At the same point the two shear moduli merge into one because of the higher symmetry of the square phase. At a still higher field we observe the reorientation transition, which occurs when $`c_r=0`$. This instability of the vortex lattice with respect to rotations also exists in the nonlocal London model mir . There are experimental indications for this instability, although its investigation has been hampered so far by the high fields required for its observation vinnikov . This reorientation transition is sensitive to the value of $`\kappa `$, and from our results we predict that this transition will not occur in the cuprates ($`\kappa 100`$), as the transition point where $`c_r=0`$ moves above $`H_{c2}`$.

We also note that anisotropic elasticity with a softer shear modulus for shearing along the sides of the square lattice than along the diagonals emerges naturally in this model. This response of the vortex lattice is measured by the ratio $`c_{66x}/c_{66,\frac{\pi }{4}}`$, which has a maximum at the transition as can be seen in Fig. 2. It has been suggested vinnikov  that this behavior explains the anisotropic orientational long range order observed in decoration patterns in LuNi2B2C, which manifests itself as a significant difference in the correlation lengths along the $`[110]`$ and $`[\overline{1}10]`$ directions.

We next turn to the effects of thermal fluctuations on the vortex lattice structure. Following Gurevich and Kogan gurevich , we will make the simplifying assumption that harmonic fluctuations renormalize the three variables $`\rho `$, $`\sigma `$, and $`\varphi `$, which characterize the structure of the vortex lattice, with the result that

$`\overline{\rho }`$ $`=\rho [1\sigma (2\rho )\overline{u^2}],`$ (14)
$`\overline{\sigma }`$ $`=\sigma [1\sigma (2\rho )\overline{u^2}],`$ (15)
$`\overline{\varphi }`$ $`=\varphi +\sigma \mathrm{cot}\varphi \overline{u^2},`$ (16)
where overbars denote thermal averages. The rhombic to square transition is then determined as the locus of points in the $`H`$$`T`$ plane where $`c_{sq}=0`$. For that purpose, the “squash” elastic modulus $`c_{sq}`$ is recalculated from Eq. (12) using the renormalized variables $`\overline{\rho }`$, $`\overline{\sigma }`$, $`\overline{\varphi }`$. This is not performed self-consistently, as we do not take into account the effects of the fluctuations on the elastic moduli themselves.

This approach is justified due to the fact that the mean squared displacement of the vortex $`\overline{u^2}`$ is finite at $`H_{\mathrm{}}(T)`$ gurevich . There are two competing parameters, the strength of the thermal fluctuations $`\chi `$ and the anisotropy parameter $`ϵ`$ which is related to the nonlocality parameter $`\rho _{GK}`$ of the nonlocal London model as $`ϵ=\left(\rho _{GK}/\xi _0\right)^2(1t^2)/12`$. To compare with the nonlocal London model results, we take gurevich  $`\lambda (T)=\lambda _0/\sqrt{1t^2}`$ and $`\xi (T)=\xi _0/\sqrt{1t^2}`$, where $`t=T/T_c`$, so that $`\chi =\chi _0t/\sqrt{1t^2}`$. The dimensionless mean squared displacement is then

$$\overline{u^2}=\frac{\chi _0\eta }{\sqrt{(1t^2)}}\frac{bt}{[(1b)^3\mathrm{ln}(1+\frac{1}{\sqrt{2b}})]^{\frac{1}{2}}},$$
(17)
with $`\eta 3`$ for LuNi2B2C gurevich . The results of the numerical solution are shown in Fig. 3, for the same parameter values used by Gurevich and Kogan gurevich , namely, $`\rho _{GK}=2.5\xi _0`$ and $`\chi _0=0.0064`$.

With fluctuations, the phase boundary curves backward and displays the “reentrant” behavior that was also obtained by Gurevich and Kogan. The physical explanation of this feature is simple: thermal fluctuations tend to wash out effects of the anisotropy, so that the preferred symmetry of the lattice is rhombic. Note that the abrupt termination of the transition line for the GL theory is an artifact of the perturbative treatment of the anisotropy; at that point the correction term $`\mathrm{\Gamma }`$ ceases to be small and the extended GL theory becomes unstable. While this shortcoming of the model could be remedied by including higher order terms in the GL free energy, the overall qualitative features of the phase diagram are expected to remain unchanged. One such choice would be an additional isotropic term. Nevertheless, the complications introduced in the analysis by the unnecessary extra parameter, would clearly outweigh any potential benefits.

We have used the GL theory with an anisotropic, symmetry breaking term to study the evolution of the vortex lattice and particularly the structural phase transition and its connection with the elastic response of the vortex lattice. We have shown that to obtain the correct phase diagram it is necessary to calculate properties up to second order in the anisotropy; at this order a rotational instability of the square vortex lattice appears. We predict that this instability should be observable in the borocarbides but not in the cuprate superconductors. Including thermal fluctuations in the model we were able to show that the rhombic to square transition line exhibits reentrant behavior, similar to what is obtained within the nonlocal London model, and in qualitative agreement with SANS studies of LuNi2B2C.

We would like to thank V. G. Kogan, A. Gurevich and P. J. Hirschfeld for useful discussions. This work was supported by NSF DMR-9978547.