text:
# 1 Introduction

## 1 Introduction

M87 is usually classified as a Fanaroff-Riley Class I (FR-I) radio galaxy having a one-sided relativistic jet  implying, in the context of unification models  that its unresolved nuclear region is a misaligned BL Lac object. The viewing angle is estimated to be $`\theta 30^{}`$ giving a maximum Doppler factor $`\delta =[\gamma _j(1\beta _j\mathrm{cos}\theta )]^1`$ of $`\delta _{\mathrm{max}}2`$ where $`\gamma _j=(1\beta _j^2)^{1/2}`$ is the jet Lorentz factor. For simplicity, we shall take $`\delta =1`$ but consider the range $`0.66\delta 1.6`$ for the M87 jet at pc scales. Although the jet emission is not strongly Doppler boosted towards us – it may even be de-boosted – its proximity to us (16.3 Mpc ) partially compensates for this, and makes TeV gamma-ray , cosmic ray  and possibly neutrino emission from this object potentially interesting.

The strong variability of the optical flux of M87 suggests that the jet emits synchrotron radiation at optical frequencies somewhere within an unresolved central region less than 5 pc in diameter . The jet has well-defined relativistic features, and it is remarkable that despite its low power the jet extends beyond 30 arcsec from the core , possibly being a remnant from when M87 was much more active. Given the large black hole mass $`3\times 10^9`$M, and the low jet luminosity $`5\times 10^{44}`$ erg s-1 , the accretion disk must be currently in a low radiative state and provide little power to the jet. As a consequence, the heating of any torus is currently inefficient and would produce little attenuation of TeV gamma-ray signals although this would change if M87 returned to an active state commensurate with its high black hole mass .

BL Lac objects, along with flat-spectrum radio quasars are collectively referred to as blazars. BL Lacs may be high-frequency or low-frequency peaked (HBLs and LBLs). Their broad-band spectra consist of two spectral components which appear as broad ‘humps’ in the spectral energy distribution (SED), and are due to emission from a jet oriented at small angle with respect to the line-of-sight. The low-energy component, is generally believed to be synchrotron emission from relativistic electrons, and extends from the radio to UV or X-ray frequencies. The origin of the high-energy component, starting at X-ray or $`\gamma `$-ray energies and extending in some cases to TeV-energies, is uncertain.

For BL Lacs which have rather weak thermal emission, the favoured “leptonic model” is the so-called synchrotron-self Compton (SSC) model in which the same relativistic electrons responsible for the low energy synchrotron hump in the SED up-scatter synchrotron photons to high energies via the Inverse Compton effect. “Hadronic models” were proposed more than 10 years ago to explain the $`\gamma `$-ray emission from blazars . Recently Mücke & Protheroe  have discussed in detail the various contributing emission processes in the synchrotron proton blazar (SPB) model. In hadronic models the relativistic jet consists of relativistic proton and electron components, which again move relativistically along the jet. High-energy radiation is produced through photomeson production, and through proton and muon synchrotron radiation, and subsequent synchrotron-pair cascading in the highly magnetized environment. In the case of BL Lacs internal photon fields (i.e. produced by synchrotron radiation from the co-accelerated electrons) serve as the target for pion photoproduction. These models can, in principle, be distinguished from leptonic models by the observation of high energy neutrinos generated in decay chains of mesons created in the photoproduction interactions (for a recent review see ).

In hadronic models, AGN would contribute also to the pool of extragalactic ultra-high energy cosmic rays (UHECR) through the decay outside the host galaxy of neutrons produced by photoproduction interactions (see  in the context of accretion shocks, and  in the context of jets) or by escape of protons directly accelerated at termination shocks of jets in giant radio lobes of Fanaroff-Riley Class II (FR-II) radio galaxies .

Greisen , and Zatsepin and Kuz’min  (GZK) showed that the nucleonic component of UHECR above $`10^{20}`$ eV will be severely attenuated in the cosmic microwave background radiation (CMBR), primarily due to pion photoproduction interactions. At $`3\times 10^{20}`$ eV the mean free path is $`5`$ Mpc and the energy-loss distance is $`20`$ Mpc  (see Stanev et al.  for recent calculations). Thus, if the UHECR are extragalactic, one would expect their spectrum to cut off at $`10^{20}`$ eV, the “GZK cut-off”. UHECR may have been observed with energies well above $`10^{20}`$ eV (see Nagano and Watson  for a review including a discussion of models for the origin of the UHECR). However, very recent data from the two largest aperture high energy cosmic ray detectors are contraditory: AGASA  observes no GZK cut-off while HiRes  observes a cut-off consistent with the expected GZK cut-off. A systematic over-estimation of energy of about 25% by AGASA or under-estimation of energy of about 25% by HiRes could account the discrepancy , and the continuation of the UHECR spectrum to energies well above $`10^{20}`$ eV is now far from certain. Nevertheless, if the spectrum does extend well beyond $`10^{20}`$ eV, even though it is only an FR-I radio galaxy M87 is an attractive possibility for the origin of the UHECR  because of its proximity.

## 2 SPB modelling M87 as a mis-aligned BL Lac

Mücke et al.  identified the critical parameters determining the properties of BL Lacs in the context of the SPB model and constructed the “average” synchrotron spectrum for each class, HBLs and LBLs, which served as the target photon distribution for their hadronic cascade. An extensive collection of blazar SEDs published by Ghisellini et al.  were used to construct the “average” SED of HBLs and LBLs by overlaying all available HBL and LBL SEDs. They found that a broken power law gave a reasonable representation of the low-frequency hump in the SED of both HBLs and LBLs

$$n(ϵ)\{\begin{array}{cc}ϵ^{\alpha _1}\hfill & \text{for}ϵ_iϵϵ_b\hfill \\ ϵ^{\alpha _2}\hfill & \text{for}ϵ_bϵϵ_c\hfill \end{array}$$

with $`\alpha _1=1.5`$, $`\alpha _2=2.25`$ and $`ϵ_i=10^{5\mathrm{}6}`$eV. The break energy $`ϵ_b`$ of LBLs varied between $``$ 0.1 eV to 1.3 eV, while the maximum synchrotron photon energy $`ϵ_c`$ ranged over two orders of magnitude, from $``$ 41 eV to 4 keV. The populated energy range of HBLs is more restricted: $`ϵ_b`$ 26 eV to 131 eV and $`ϵ_c`$ 4.1 keV to 41 keV. The peak of the low-energy SED was $`\mathrm{log}\nu L_\nu ^{\mathrm{max}}(\mathrm{erg}/\mathrm{s})45.646.1`$ for LBLs and $`\mathrm{log}\nu L_\nu ^{\mathrm{max}}(\mathrm{erg}/\mathrm{s})43.443.8`$ for HBLs. Mücke et al.  defined the “average LBL” by $`ϵ_b=`$1.3 eV, $`ϵ_c=`$4.1 keV and $`\mathrm{log}\nu L_\nu ^{\mathrm{max}}\text{ (erg/s)}=46.1`$, and the “average HBL” by $`ϵ_b=`$131 eV, $`ϵ_c=`$41 keV and $`\mathrm{log}\nu L_\nu ^{\mathrm{max}}=43.8`$. The parametrizations for HBLs and LBLs are visualized in Fig. 1, with the dashed lines showing the “average” SED and the hatched regions showing the range of the SEDs for the two cases.

Assuming that the unresolved nuclear emission in M87, i.e. from the pc scale jet, has $`\delta =1`$ we have added to Fig. 1 the radio , infrared  and X-ray  data on the nuclear region of M87 Doppler boosted by $`\delta =10`$ to mimic how M87 would appear if its jet were closely aligned towards us such that it would be classed as a BL Lac rather than a FR-I. Since the Doppler factor of M87 is not well known, taking the range from 0.66 to 1.5, the observed SED as viewed along the axis with $`\delta =10`$ would move along the error bars added to the data. From Fig. 1, the rather high radio-infrared luminosity would suggest that M87 could be a mis-aligned LBL, in which case the X-ray emission would be part of the high-energy hump in the SED due to cascading initiated by photons and electrons from pion decay or by synchrotron photons emitted by protons and muons. However, the X-ray emission is at a level closer to that expected from an HBL, in which case it would be mainly due to synchrotron emission by directly accelerated electrons. One should bear in mind that there is considerable variation in the SEDs of HBLs and LBLs, and so one should not expect the SED of an individual BL Lac to be identical to our “average” spectrum. This is particularly true for the M87 nuclear region where the observations at radio, infrared and X-rays were not simultaneous and are known to vary. In addition, Wilson and Yang  have calculated from the Chandra X-ray flux and spectrum observed from the inner part of the jet that there is intrinsic absorption by cold matter with an equivalent hydrogen column density of $`3`$$`5\times 10^{20}`$ cm-2 in M87, and so the unattenuated X-ray flux could be higher than plotted. Thus we shall be satisfied if the observed SED is within a factor $`2`$ of an “average” SED.

In Fig. 2 we show the SEDs of HBLs modelled by Mücke et al.  in the context of the SPB model as they would appear if observed at $`\delta =1`$ as FR-I radio galaxies at a distance of 16 Mpc appropriate to M87. Results taken from Mücke et al.  for four peak luminosities of the low energy hump of the SED are shown after shifting to $`\delta =1`$. As we see, there is little variation in the high energy hump which peaks at $`10`$$`100`$ GeV at a level just below the sensitivities of EGRET and the Whipple High Energy Gamma Ray Telescope, but well above the sensitivity of GLAST and future northern hemisphere large-area atmospheric Cherenkov telescopes (ACT) such as MAGIC and VERITAS (sensitivities for all these telescopes taken from ref.  are shown together with the upper limit from Whipple ). One should bear in mind, however, that the uncertainty in the Doppler factor of M87 will give rise to a much larger uncertainty in the bolometric luminosity. These uncertainties are indicated in the theoretical SED and neutrino spectra by the slanted error-bars. We have extended the work of Mücke et al.  by calculating the neutrino output for these four cases and plot the expected flux of muon-neutrinos ($`\nu _\mu `$+$`\overline{\nu }_\mu `$). The highest luminosity modelled for the low-energy hump of the SED (solid curve) is in closest agreement with the M87 SED and predicts the highest neutrino flux. We have added an upper-limit for the AMANDA-B10 neutrino telescope by taking the diffuse neutrino background limit  and multiplying by $`2\pi `$ sr (neutrino telescopes detect upward-going neutrinos). Sensitivities calculated by Albuquerque et al.  of AMANDA-II and IceCube to $`E^2`$ diffuse intensity, again multiplied by $`2\pi `$ sr to convert to point source sensitivity have also been added. Unlike the 100 GeV gamma-ray emission which should easily be detected by future telescopes, the neutrino flux is well below detection levels of future neutrino telescopes.

In Fig. 3 we model M87 as a mis-aligned LBL in the context of the SPB model. In this case, we see that the gamma-ray emission cuts off at 10 GeV–1 TeV depending on the luminosity of the low energy hump of the SED. This is because for the higher target photon densities in LBL, pion photoproduction losses determine the maximum proton energy, and hence the maximum gamma-ray energy which results from the cascades initiated by proton synchrotron radiation and pion decay (including pion and muon synchrotron radiation). This also affects the maximum neutrino energy in the same way. The lowest luminosity modelled for the low-energy hump of the SED (chain curve) is in closest agreement with the M87 SED and predicts the highest gamma-ray and neutrino energies. The gamma ray flux is a factor $`3`$–10 lower than for the HBL case, but is still detectable by future telescopes. As with the HBL case, the neutrino flux is below the sensitivity of IceCube. However, it is not impossible that in the case of rapid flaring, M87 might just be detected as it is at a declination of $`+`$$`12^{}`$ and is ideally located for observation by a giant neutrino telescope to be located at the South Pole such as IceCube because it is sufficiently below the horizon to eliminate cosmic ray events, while not being so far below the horizon that neutrino absorption by the Earth’s core becomes important.

## 3 M87 as a source of UHECR

Energetic protons magnetically trapped in AGN jets lose energy predominantly by Bethe-Heitler pair production and pion photoproduction on ambient radiation fields, or by adiabatic deceleration as the jet expands. Neutron production in pion photoproduction sources (e.g. $`p\gamma n\pi ^+`$) provides a mechanism for escape of cosmic rays from the jet. Neutrons decay typically after travelling $`(E_n/10^{20}`$eV) Mpc, which for UHECR is well outside the host galaxy. Recently Ahn et al.  and Biermann et al.  have shown that in the assumption of a Parker-spiral galactic wind magnetic structure out to $`1.5`$ Mpc, the galactic wind of our Galaxy poses no restriction to the entry to our Galaxy of UHECR from the general direction of M87, and may even lead to clustering in the arrival directions of cosmic rays from M87 observed from Earth. These preferred arrival directions and clustering may, however, be mainly due to the Galactic wind topology and contain little information about the direction of cosmic ray sources .

The observed intensity of UHECR  (multipled by energy squared) is plotted in Fig. 4. For the most promising LBL model (chain curve in Fig. 3) we plot (chain curve) the flux of neutrons divided by $`4\pi `$ sr (to convert flux to average intensity) that would be observed at Earth from M87 if the neutrons did not decay (which of course they do). Since these neutrons would travel to Earth in straight lines, and since we take M87 to have $`\delta (\theta 30^{})=1`$, the maximum neutron energy that would arrive at Earth would be approximately the maximum jet-frame energy of the accelerated protons, which would be $`E_n=3\times 10^{19}`$ eV for the model shown. However, neutrons decay into protons whose directions are isotropised in the intergalactic magnetic field before and during propagation to Earth. Hence we would have contributions to the protons arriving at Earth from neutrons emitted at all angles with respect to the jet axis, and hence from the full range of Doppler factors corresponding to the jet Lorentz factor, i.e. approximately $`1/2\gamma _j2\gamma _j`$. For example, if the jet Lorentz factor were $`\gamma _j=5`$ a neutron with jet-frame energy $`E_n^{}=3\times 10^{19}`$ eV emitted along the jet axis ($`\theta =0^{}`$) would decay to a proton with galaxy-frame energy $`E3\times 10^{20}`$ eV. What we would observe therefore depends on the galaxy-frame angle-averaged neutron luminosity on emission, and would also have contributions from both jets (i.e. twice that for one jet). Given that the chain curve in Fig. 4 is for $`\delta (\theta 30^{})=1`$, it is related directly to the jet-frame bolometric luminosity $`E_{n}^{}{}_{}{}^{2}\dot{N}_n^{}(E_n^{})`$. For electromagnetic radiation and neutrinos the observed bolometric luminosity is $`\nu L_\nu ^{\mathrm{obs}}(\theta )=\delta ^4(\theta )\nu ^{}L_\nu ^{}^{}`$ and the angle-averaged emitted bolometric luminosity is $`\nu L_\nu ^{\mathrm{em}}=(4\pi )^1\delta ^3(\theta )\gamma _j^1\nu ^{}L_\nu ^{}^{}𝑑\mathrm{\Omega }`$. Hence, by analogy, for ultra-relativistic protons of energy $`EE_n`$ from neutron decay, and remembering that both jets contribute,

$$E^2\dot{N}(E)=\frac{2}{4\pi }\delta ^3(\theta )\gamma _j^1E_{n}^{}{}_{}{}^{2}\dot{N}_n^{}(E_n^{})𝑑\mathrm{\Omega }.$$
(1)
We have added to Fig. 4 (solid curve) the resulting cosmic ray flux divided by $`4\pi `$ sr assuming $`\gamma _j=5`$ and straight-line propagation, i.e.

$$E^2I(E)=\frac{E^2\dot{N}(E)}{(4\pi )4\pi d^2}$$
(2)
corresponding to number density

$$n(E)=\frac{\dot{N}(E)}{4\pi d^2c}.$$
(3)
At $`5\times 10^{19}`$$`3\times 10^{20}`$ eV this is a factor $`20`$ below the observed UHECR.

Of course, so far this neglects effects of diffusive propagation to Earth and interactions with the CMBR. We shall consider first the effects of diffusion. For the case of a constant injection rate of cosmic rays, and an infinite homogeneous diffusing medium with diffusion coefficient $`D(E)`$ the number density at distance $`d`$ is simply

$$n(E)=\frac{\dot{N}(E)}{4\pi D(E)d}.$$
(4)
If we take $`D(E)=\lambda (E)c/3`$ then the observed intensity will be enhanced with respect to straight-line propagation by a factor

$$g(E)=dc/D(E)=3d/\lambda (E).$$
(5)
For a simple model of intergalactic space with a low magnetic field ($`10^9`$ G) and a cell size $`L_C>100`$ kpc, taking $`\lambda (E)L_C`$ would give an enhancement factor up to $`g(E)<500`$, such that the predicted UHECR from M87 could explain the observed intensity above $`10^{19}`$ eV. While diffusion appears to increase the cosmic ray intensity from M87, it also increases the travel time or effective distance travelled from M87 to Earth,

$$d_{\mathrm{eff}}=d^2c/2D(E)=dg(E)/2.$$
(6)
Hence for an enhancement factor of $`g=20`$ the effective distance is $`d_{\mathrm{eff}}=160`$ Mpc. At $`10^{19}`$ eV the mean energy-loss distance is approximately 1 Gpc and at $`10^{20}`$ eV the mean energy-loss distance is approximately 160 Mpc, and so in this case there would only be a significant energy loss above $`10^{20}`$ eV.

The situation could actually be quite different if we consider an intergalactic magnetic field (IGMF) structure based on the observation of microgauss fields in clusters of galaxies , and of clusters occurring in networks of “walls” separated by “voids”. It would then be reasonable to expect relatively high fields of $`10^7`$$`10^6`$ G in the walls, and much lower fields, $`10^{11}`$$`10^9`$ G, in the voids. Recently, Stanev et al.  have considered propagation of UHECR in three different models of the IGMF in the local supergalactic structure in which the regular field in the high-field region was $`10^8`$ G. We consider propagation in a simple wall/void model similar to that used by Medina Tanco  to illustrate how UHECR would be deflected in complex IGMF structures. As M87 is at the centre of the Virgo Cluster and is close to the super-galactic plane we assume it is embedded in the higher field region. To illustrate the effect of propagation from M87 in a wall/void type IGMF we define the origin of coordinates to be at M87 with the mid-plane of the wall corresponding to the $`x`$$`y`$ plane, the wall occupying $`|z|<2.5`$ Mpc and the void occupying $`|z|>7.5`$ Mpc. Following Medina Tanco  we adopt a regular magnetic field of $`10^7`$ G in the $`x`$-direction in the wall and $`10^{10}`$ G in the void, with a transition region sandwiched between the wall and the void in which the magnetic field drops exponentially from $`10^7`$ G to $`10^{10}`$ G. The magnetic field contains an irregular component having $`|\stackrel{}{B}_{\mathrm{irreg}}|^2^{1/2}`$ equal to 30% of the regular component, and has a Kolmogorov spectrum of turbulence. The irregular field is modelled as described in ref.  with wavenumbers corresponding to turbulence scales $`L_C/2^n`$ with $`L_C=2.5`$ Mpc and $`n=0,\mathrm{},3`$. Particles are injected isotropically at the origin and energy losses due to Bethe-Heitler pair production and pion photoproduction are included. A typical simulation for $`E_0=10^{20}`$ eV is shown in Fig. 5.

Whereas Medina Tanco  was most interested in the angular deflection of arriving cosmic rays relative to the direction of the cosmic ray source, we are more interested in the travel times of cosmic rays to Earth from the distance $`d16`$ Mpc to M87, and the total amount of time spent within unit volume at distance $`d`$ from M87. We divide the spherical shell corresponding to 15 Mpc $`<d<17`$ Mpc into 2160 equal cells of volume $`V`$ with a grid in spherical coordinates $`(\theta ,\varphi )`$ using 72 equally spaced $`\varphi `$-values and 30 equally spaced $`\mathrm{cos}\theta `$-values. Particle orbits are advanced in time steps of size $`\mathrm{\Delta }t=0.03`$ Mpc/$`c`$ chosen to be much smaller than the smallest magnetic structure simulated. At each step, the proton energy is reduced according to the average energy-loss rate (see e.g. figure 1 of ref. ). A check is made at the end of each time step to determine whether the proton is within the spherical shell, and if it is then $`\mathrm{\Delta }t`$ is added to the total time spent in the cell corresponding to its spherical coordinates $`(\theta ,\varphi )`$ – i.e. if on $`k`$ occasions a particle was found to be located inside the $`i`$th cell at the end of a time step, then it would have spent approximately a total time $`k\mathrm{\Delta }t`$ inside the $`i`$th cell since being injected at the origin. Dividing by the number of simulations one obtains the average time spent in a particular cell per injected proton $`\tau (\theta ,\varphi )`$. Similarly, the time spent travelling since injection at M87 is added to the ages of the protons “detected” in that cell. Dividing by the total number of occasions a particle was “detected” in that cell gives the average age $`\overline{t}(\theta ,\varphi )`$ of the protons present in that cell.

Results based on the average time per unit volume in $`10^6`$ simulations are shown for two initial energies $`10^{19}`$ eV and $`10^{20}`$ eV in Fig. 6. The gray-scale shows the enhancement factor

$$g(\theta ,\varphi )=4\pi d^2c\tau (\theta ,\varphi )/V.$$
(7)
The strong peaks at $`(0^{},0)`$ and $`(180^{},0)`$ show that despite the presence of the turbulent field component, cross-field diffusion is not strong enough to give rise to significant fluxes far away from the regular field threading the source, at least for source distances as small as 16 Mpc for this field topology, and we note that Stanev et al.  arrived at a similar conclusion about cross-field diffusion. The peaks in $`g(\theta ,\varphi )`$ shown in Fig. 6 are at a level $`g(\theta ,\varphi )10^3`$ over a disk of radius $`5^{}`$ ($`1.4`$ Mpc) outside of which $`g(\theta ,\varphi )1`$, indicating that if the magnetic field topology is such that the field lines connect the vicinity of M87 (within $`1`$ Mpc) to our Galaxy (within $`1`$ Mpc) then cosmic rays from M87 should clearly be observed. The average ages are $`\overline{t}200`$ Mpc/c ($`E_0=10^{19}`$eV) and $`\overline{t}100`$ Mpc/c ($`E_0=10^{20}`$eV) indicating their final energies would be $`8\times 10^{18}`$eV and $`5\times 10^{19}`$eV respectively.

## 4 Conclusion

We find the unresolved nuclear core of M87 to be consistent with it being a mis-aligned BL Lac object. In the context of the Synchrotron Proton Blazar model, M87 could be either an HBL or an LBL. In both cases, we predict gamma-ray emission at levels detectable by GLAST and MAGIC, and possibly by VERITAS. Neutrino detection is not expected, except possibly during an extreme flare in the LBL case.

Ultra-high energy neutrons produced in pion photoproduction interactions will escape from the host galaxy where they decay into protons. Because cosmic rays are deflected in the IGMF, protons resulting from neutrons emitted in all directions can contribute to the observed UHECR. Hence, even though the electromagnetic radiation we observe from the M87 jet is not significantly Doppler boosted in energy, the cosmic ray output will be. By this mechanism, the unresolved nuclear core of the M87 jet could emit UHECR with energies up to at least $`3\times 10^{20}`$ eV.

We predict the UHECR output from M87 to be at a level such that if UHECRs travelled in straight lines they would give an average intensity at Earth a factor $``$20 below that observed. For a constant UHECR output and simple isotropic diffusion with a scattering mean free path $`\lambda `$ less than the distance $`d`$ to M87, the predicted intensity would be enhanced by a factor $`3d/\lambda `$ relative to straight line propagation. This enhancement is accompanied by increased travel times, and so higher energy losses in the CMBR. Nevertheless, taking reasonable scattering mean free paths we find that for simple diffusion models protons from decay of ultra-high energy neutrons produced by pion photoproduction in the M87 jet could easily account for all the observed UHECR.

In perhaps more realistic wall/void models of IGMF structure, M87 would only be a source of the observed UHECR if the topology of the IGMF between M87 and our Galaxy is favourable. Note, however, that because of its very high black hole mass M87 was probably much more active at earlier times than at present. Many objects exhibit a high state for $``$5% of the time, and since we estimate the average travel time of UHECR from M87 to Earth to be a factor $``$5–10 times the light propagation time it is indeed possible that the UHECR observed now were emitted when M87 was in a high state and possibly a more powerful cosmic ray source. Similarly, if M87 was in the recent past a powerful FR-II radio galaxy, it could have been possible to accelerate UHECR at shocks in the giant radio lobes where the jets terminate as in the model of Rachen and Biermann , and these could also contribute to the observed UHECR.

## Acknowledgments

We thank Todor Stanev for very helpful comments. This work was supported by a Discovery Project grant from the Australian Research Council and a grant from the University of Adelaide to RJP. AR thanks the Bundesministerium für Bildung und Forschung for financial support through DESY grant Verbundforschung 05CH1PCA/6.


format: text

text:
# Vortex-Line Phase Diagram for Anisotropic Superconductors

## I Introduction

One important theoretical problem in the study of type-II superconductors is the zero temperature equilibrium vortex phase diagram for clean, uniaxially anisotropic, bulk materials subjected to an external magnetic field, $`𝐇`$ tilted with respect to the axis of anisotropy ($`c`$-axis). In the case of isotropic materials the vortex phase is the well known triangular vortex-line lattice with the vortex lines oriented parallel to $`𝐇`$. In anisotropic materials the vortex phase is still a vortex-line lattice but, in general, the lattice is not triangular nor are the vortex lines oriented along $`𝐇`$ .

Previous work find that the vortex phases for the anisotropic superconductor are related to those for the isotropic one when $`𝐇`$ is oriented along the c-axis, in which case they are identical, when $`𝐇`$ is along the a-b plane or for $`H`$ large compared with the lower critical field . In the latter two cases, the vortex-line-lattices are related to the triangular lattice by scaling and the vortex lines are essentially parallel to $`𝐇`$. For general orientations of $`𝐇`$ this is not the case. The vortex phases were first obtained by Buzdin and Simonov  by numerically minimizing the Gibbs free-energy in the London limit with material parameters appropriate for YBCO. Our calculation, like their’s, minimizes the London Gibbs free-energy numerically. However, we go beyond the work of these authors by incorporating into our calculations recent results for the vortex state at the lower critical field, not taken into account by them . We show that this leads to new vortex phases close to the lower critical field.

It was found first by Grishin, Martynovich and Yampol’skii  and later by others , that a pair of straight vortex lines parallel to each other and coplanar with respect to the c-axis attract at large distances. A consequence of this is that a chain of such lines has lower energy than the same vortex lines placed far apart. On the basis of this result Buzdin and co-workers  suggested that at the lower critical field the mixed state consists of a vanishing small density of vortex line chains instead of vortex lines. However this argument does not apply in general. As pointed out in Ref., in order to determine whether vortex lines or vortex-line chains constitute the vortex phase at the lower critical field it is necessary to minimize the Gibbs free-energy for each one of them and obtain their respective lower critical fields. The one with the lower critical field is the true equilibrium phase. One important aspect of the calculation of Ref. is the use of an expression for the tilted vortex-line self energy with an appropriate anisotropic core cutoff function. This expression was derived by Subdo and Brandt  and later used in Ref. to show, by minimizing the London Gibbs free-energy for a vanishing small density of vortex lines, that for anisotropy parameters typical of BSCCO there is coexistence two such states, with different vortex lines tilt angles, for one particular orientation of $`𝐇`$. These authors also showed that if a self-energy expression based on an isotropic core cutoff, such as the one adopted by Buzdin and Simonov , is used instead, the coexistence disappears. It is found in Ref. that at the lower critical field and for anisotropy parameters typical of BSCCO the equilibrium vortex state for $`𝐇`$ tilted by $`\alpha <\alpha _c=7.9^o`$ with respect to the c-axis is a dilute vortex-line chain state whereas for $`\alpha >\alpha _c=7.9^o`$ it is a dilute vortex-line one. For anisotropy parameters typical of YBCO Ref. finds that the lower critical fields for vortex lines and vortex-line chains differ so little from each other that it is not possible do decide which one corresponds to true equilibrium.

In this paper, we minimize the exact Gibbs free-energy in the London limit using Sudbo and Brandt  expression the vortex line self energy. We assume that the equilibrium state the vortex lines are straight and parallel to each other and the VLL has only one vortex line per unit cell, and obtain the vortex lines orientation and VLL shape as functions of $`𝐇`$, by minimizing the Gibbs free-energy in the London limit.

The free parameters in our calculation for given $`𝐇`$, specified here by its magnitude, $`H`$, and its tilt angle with respect to the c-axis, $`\alpha `$, are the vortex lines tilt angle, $`\theta `$, and the VLL primitive unit cell vectors $`𝐋_\mathrm{𝟏}`$ and $`𝐋_\mathrm{𝟐}`$, defined in the plane perpendicular to the vortex lines. Since the vortex lines are straight and parallel to each other, the magnetic induction $`𝐁`$ is also tilted by $`\theta `$ with respect to the c-axis and, as explained in detail in Sec. II A, it is coplanar with $`𝐇`$ and the c-axis. The orientations of $`𝐇`$, $`𝐁`$, and the c-axis, and the VLL primitive unit cell are shown in Fig. 1. Hereafter we specify the VLL unit cell by the magnitudes of $`𝐋_\mathrm{𝟏}`$ and $`𝐋_\mathrm{𝟐}`$, denoted by $`L_1`$ and $`L_2`$, by the angle between $`𝐋_\mathrm{𝟏}`$ and $`𝐋_\mathrm{𝟐}`$, $`\phi `$, and by and the angle of rotation of the VLL unit cell with respect to $`𝐁`$, $`\psi `$, as shown in Fig. 1.

As discussed in Sec. II A, a minimum of the Gibbs free-energy occurs for $`\psi =0`$, where $`\psi `$ is defined in Fig. 1. We assume that this is the absolute minimum and carry out the calculations with $`\psi =0`$. This means that $`𝐋_\mathrm{𝟏}`$ is also coplanar with $`𝐇`$ and the c-axis,oriented along the x-direction defined in Fig. 1. To simplify the calculations we set $`\phi =90^o`$. This means that the VLL unit cell is a rectangle, with $`𝐋_\mathrm{𝟐}`$ perpendicular to the plane defined by $`𝐇`$, $`𝐁`$ and the c-axis (Fig. 1). We find that this simplifications does not change our main results. Under these conditions we have only three free parameters, to be determined by minimizing the Gibbs free-energy: $`\theta `$, $`L_1`$ and $`L_2`$. This minimization is carried out numerically, using a fast convergent series to calculate lattice sums entering the expression for the VLL interaction energy . This series, together with the small number of free parameters, allow us to create an efficient and accurate numerical algorithm.

We study in detail the phase diagram for anisotropy parameters typical of YBCO and BSCCO. For YBCCO our results are essentially the same as those obtained by Buzdin and Simonov, except very close to the lower critical field. For BSCCO our results are new.

The main results of this paper are summarized in the generic phase diagram show in Fig. 2. We find three distinct phases, hereafter called phase-I, phase-II, and phase-III. We find that phase-I is new and closely related to the vortex-line equilibrium state at the lower critical field. We also find that phase- II is essentially identical to that found by Buzdin and Simonov and that phase-III is that expected for fields much larger than the lower critical field. In this limit the VLL unit cell is related to that appropriate for the isotropic superconductor by scaling, and $`𝐁𝐇`$ ($`\theta \alpha `$.

The lower critical field in isotropic superconductors is independent of the orientation $`(\alpha )`$, and has magnitude $`H_{c1}`$. For $`H`$ just above $`H_{c1}`$, the vortex state corresponds to vortex lines parallel to $`𝐇`$ very far apart. In our notation this state is described by $`\theta =\alpha `$ and $`L_1`$, $`L_2=\mathrm{}`$. In a superconductor with uniaxial anisotropy, the magnitude of the lower critical field depends on $`\alpha `$, hereafter denoted by $`H_{c1}(\alpha )`$ . The vortex state for $`H`$ just above $`H_{c1}(\alpha )`$ consists of vortex lines not parallel to $`𝐇`$ ($`\theta \alpha `$) in one of two possible configurations, named dilute vortex line state (DVL) and dilute vortex-line chain state (DVLC) . The DVL is similar to the traditional vortex state for the isotropic superconductor, in the sense that $`L_1`$, $`L_2=\mathrm{}`$, whereas the DVLC corresponds to vortex-line chains along the x-direction with period $`L_1`$, very far apart, that is $`L_1`$ finite and $`L_2=\mathrm{}`$. The latter possibility results from attractive interactions between vortex lines that only exist in anisotropic superconductors . To decide if the DVL or the DVLC is the equilibrium vortex state it is necessary to calculate $`H_{c1}(\alpha )`$ for each one. The equilibrium configuration is that with the smallest $`H_{c1}(\alpha )`$. This calculation is reported in Ref.. The results for anisotropy parameters typical of BSCCO and YBCO, referred to as strong an moderate anisotropies, respectively, are shown in Fig. 3. It is found that for both anisotropies there is a value of the external field tilt angle, $`\alpha _c`$, around which the equilibrium vortex states undergoes a significant change. For $`\alpha >\alpha _c`$ and both for strong and moderate anisotropies the equilibrium state is a DVL with vortex lines tilted at $`\theta 90^o`$, that is parallel to the ab-plane. For strong anisotropy the vortex state for $`\alpha <\alpha _c=7.9^o`$ is a DVLC with vortex lines tilted at $`\theta <60^o`$. At $`\alpha =\alpha _c`$, there is coexistence of a DVL and DVLC. A discontinuous change in $`\theta `$ takes place at $`\alpha =\alpha _c`$. For a moderate anisotropy, $`\theta `$ is found to change smoothly with $`\alpha `$ and $`\alpha _c20^o`$. The calculation on Ref. finds that for moderate anisotropy the DVL and DVLC values of $`H_{c1}(\alpha )`$ are so close that their difference is within the numerical errors of the calculation.

We find that phases I and II evolve, respectively, from the $`\alpha <\alpha _c`$ and $`\alpha >\alpha _c`$ vortex states at $`H=H_{c1}(\alpha )`$, represented by the horizontal line in Fig. 2.

Phase-II is essentially an extension into the $`H>H_{c1}(\alpha )`$ region of the phase diagram of the vortex state found for $`\alpha >\alpha _c`$ at $`H=H_{c1}(\alpha )`$. Its predominant characteristic is that $`\theta 90^o`$ and, consequently, the VLL within phase-II is essentially identical to that for vortex lines parallel to the ab-plane. In this case the unit cell is, approximately, related to that for the isotropic superconductor by scaling relations that depend only on the anisotropy. For strong and moderate anisotropies the unit cell is a rectangle elongated in the y-direction, that is $`L_1L_2`$. Phase-II is limited from above by the line $`H=H^{}(\alpha )`$, which gives the field value around which $`\theta `$ changes from $`90^o`$ to a lower value, accompanied by significant changes in the VLL unit cell. At $`\alpha =\alpha _c`$, $`H^{}(\alpha )=H_{c1}(\alpha )`$.

Phase-I is the region of the phase diagram where the VLL properties are strongly influenced by the attractive interactions between pairs of vortex lines aligned along the x-directions. For both strong and moderate anisotropies we find that phase-I is not limited to $`\alpha <\alpha _c`$, but it also exists in a small part of the region $`\alpha >\alpha _c`$ and $`H>H^{}(\alpha )`$. For strong anisotropy we find that phase-I consists of vortex-line chains along the x-direction, with a period $`L_1`$ that is nearly identical to that of an isolated chain with the same tilt angle (which differs from that of an isolated vortex), separated from each other by $`L_2L_1`$. This means that the attractive interactions between vortex lines prevail in determining the VLL periodicity along the x-direction. We also find that as $`HH_{c1}(\alpha )`$ from above, the equilibrium VLL approaches the same DVLC found in Ref. for $`H=H_{c1}(\alpha )`$. For moderate anisotropy we obtain a new result for the vortex state at $`H=H_{c1}(\alpha )`$. By extrapolating our data as $`HH^{}(\alpha )`$ from above, we conclude that the vortex state at $`H=H_{c1}(\alpha )`$ is a DVLC. However, within phase-I for moderate anisotropy, the attractive interactions between vortex lines are so weak that vortex-line chains along the x-direction with periods close to those of isolated chains only exist within phase-I near $`H=H_{c1}(\alpha )`$. However, distinctive behavior, due to attractive vortex-line interactions, can still be identified within phase-I.

This paper is organized as follows. In Sec. II A we review London theory for the energy of the vortex-line system under consideration and the known results for the phase diagram at $`H=H_{c1}(\alpha )`$ and for particular cases at $`H>H_{c1}(\alpha )`$. In Sec. II B we discuss the numerical method used to minimize the Gibbs free energy. In Sec. III we report the results of this minimization, the detailed properties of each phase, and discuss how the phases are identified. Our conclusions are stated in Sec. IV. In the Appendix the relevant formulas for the fast convergent series used in our in our calculation are summarized.

## II Method of Calculation

Here, we review London theory for this problem, its known results for the phase diagram, and describe some details of our numerical method.

### A London Theory

To obtain the zero-temperature phase diagram for a superconductor in an applied field it is necessary to determine the vortex arrangement that minimizes the Gibbs free-energy density

$$G=F\frac{𝐇𝐁}{4\pi }$$
(1)
with $`𝐇`$ and the volume kept constant. In Eq. (1) $`F`$ is the vortex energy density and $`𝐁`$ is the magnetic induction. We assume that the zero-temperature equilibrium vortex configuration consist of straight vortex lines, parallel to each other, arranged on a periodic VLL, with a single vortex line per primitive unit cell. In this case $`𝐁`$ is parallel to the vortex lines direction and has magnitude $`B=\mathrm{\Phi }_0/A_c`$, where $`A_c`$ is the VLL primitive unit cell area. Minimization of $`G`$ gives then the vortex lines orientation and unit cell shape that corresponds to equilibrium. Because $`F`$ is unchanged, if the vortex lines are rotated around the $`c`$-axis, the equilibrium $`𝐁`$ lies in the same plane as $`𝐇`$ and the $`c`$-axis. Thus, for $`𝐇`$ with magnitude $`H`$ and tilted with respect to the $`c`$-axis by an angle $`\alpha `$, $`G`$ depends on the vortex lines tilt angle with respect to the $`c`$-axis, $`\theta `$, and on the parameters defining the VLL primitive unit cell. In the system of axis defined in Fig. 1.a, the generic VLL primitive unit cell lies in the $`xy`$ plane, and is defined by the primitive unit vectors $`𝐋_\mathrm{𝟏}`$ and $`𝐋_\mathrm{𝟐}`$, as shown in Fig. 1.b. In this case $`A_c=L_1L_2\mathrm{sin}\phi `$.

In the London limit the superconductor is characterized by the penetration depths $`\lambda _{ab}`$ and $`\lambda _c`$, for currents parallel to the $`ab`$-plane and to the $`c`$-direction, respectively. The free energy density of a generic arrangement of these vortex lines can be written as the sum of pairwise interactions

$$F=\frac{\mathrm{\Phi }_0^2}{8\pi A}\underset{i,j}{}f(𝐫_i𝐫_j),$$
(2)
where $`𝐫_i`$ is the i-th vortex-line position vector in $`xy`$ plane, $`A`$ is the sample area in this plane, and $`f(𝐫)`$ is the Fourier transform of

$$f(𝐤)=e^{2g(𝐤)}\frac{1+\lambda _\theta ^2k^2}{(1+\lambda _{ab}^2k^2)(1+\lambda _\theta ^2k_x^2+\lambda _c^2k_y^2)},$$
(3)
where

$$\lambda _\theta ^2=\lambda _{ab}^2\mathrm{sin}^2\theta +\lambda _c^2\mathrm{cos}^2\theta ,$$
(4)
and $`g(𝐤)`$ is the vortex core cutoff function . It is costumery to write $`G`$ as

$$G(\theta ,𝐋_\mathrm{𝟏},𝐋_\mathrm{𝟐};𝐇)=\frac{B}{\mathrm{\Phi }_0}[\epsilon _{\mathrm{sf}}(\theta )+\epsilon _{\mathrm{in}}(\theta ,𝐋_\mathrm{𝟏},𝐋_\mathrm{𝟐})\frac{\mathrm{\Phi }_0}{4\pi }H\mathrm{cos}(\theta \alpha )],$$
(5)
where $`\epsilon _{\mathrm{sf}}(\theta )`$ is the vortex line self-energy and $`\epsilon _{\mathrm{in}}(\theta ,𝐋_\mathrm{𝟏},𝐋_\mathrm{𝟐})`$ is the interaction energy per vortex line. In terms of $`f(𝐫)`$, defined in Eq. (2), these quantities are given by

$$\epsilon _{sf}(\theta )=\frac{\mathrm{\Phi }_0}{8\pi }f(𝐫=0),$$
(6)
and

$$\epsilon _{in}=\frac{\mathrm{\Phi }_0}{8\pi }\underset{ij}{}f(𝐑_i𝐑_j),$$
(7)
where $`𝐑_i`$ denotes the VLL positions ($`𝐑_i=n_1𝐋_1+n_2𝐋_2`$, $`n_1`$, $`n_2=`$ integer).

The minima of $`G`$, Eq. (5), are known in the following cases.

i) Lower critical field: $`H=H_{c1}(\alpha )`$. The equilibrium vortex phase at the lower critical field is studied by Sudbø, Brandt and Huse , assuming that it consists of a dilute arrangement of straight vortex lines (DVL). Accordingly, these authors minimize $`G`$, Eq. (5), neglecting $`\epsilon _{\mathrm{in}}`$, and using for $`\epsilon _{\mathrm{sf}}`$ the expression derived by Sudbø, and Brandt , based on an elliptic core cutoff function,

$$\epsilon _{\mathrm{sf}}(\theta )=\epsilon _0\frac{\lambda _\theta }{\lambda _c}[\mathrm{ln}(\kappa /\gamma )+\frac{\lambda _c^2\mathrm{cos}^2\theta }{\lambda _c^2\mathrm{cos}^2\theta +\lambda _\theta ^2}\mathrm{ln}\frac{\gamma ^2(\lambda _c^2+\lambda _\theta ^2)}{2\lambda _\theta ^2}],$$
(8)
where $`\epsilon _0=(\mathrm{\Phi }_0/4\pi \lambda _{ab})^2`$, $`\gamma =\lambda _{ab}/\lambda _c=\xi _c/\xi _{ab}`$, $`\kappa =\lambda _{ab}/\xi _{ab}`$, $`\xi _{ab}`$ being the $`ab`$-plane coherence length. This line energy is obtained from the anisotropic Ginzburg-Landau theory using the Klemm-Clem transformation . The elliptical core cutoff has semi major axis $`\xi _{ab}^1`$ and semi minor axis $`\xi _\theta ^1`$ with $`\xi _\theta ^2=\xi _{ab}^2\mathrm{cos}^2\theta +\xi _c^2\mathrm{sin}^2\theta `$. Sudbø, Brandt and Huse  obtain $`H_{c1}(\alpha )`$ for several values of the anisotropy parameter $`\gamma `$ and of $`\kappa `$, and find that for certain ranges of $`\gamma `$ and $`\kappa `$, there is a value of $`\alpha `$, $`\alpha =\alpha _c`$, for which coexistence of two DVL differing from each other by the equilibrium $`\theta `$ takes place.

In a recent publication by some of us , the Sudbø, Brandt and Huse  calculation was generalized to account for two possible vortex phases at the lower critical field: a dilute vortex-line state (DVL) considered by them and a dilute vortex-line chains arrangement (DVLC). The DVLC consists of vortex lines tilted with respect to the $`c`$-axis, parallel to one another, and aligned along the x-direction, forming a periodic chain. This is possible because the interaction between a pair of such vortices is attractive at large distances . In Ref:, $`G`$, Eq. (5), is minimized in the limit of vanishing vortex density, for each one of these possibilities. For the DVL, $`\epsilon _{\mathrm{int}}`$ is neglected in Eq. (5), whereas for the DVLC $`\epsilon _{\mathrm{int}}`$ in Eq. (5) is identified with the interaction energy per vortex line of an isolated chain. The lower critical field as a function of the tilt angle $`\alpha `$, $`H_{c1}(\alpha )`$ is then obtained for each possibility, and the equilibrium vortex phase at $`H=H_{c1}(\alpha )`$ is identified as that with the smallest $`H_{c1}(\alpha )`$.

In Ref:, two sets of anisotropy parameters are considered:

Strong anisotropy: $`\kappa =10`$, $`\gamma =1/\sqrt{200}`$

Moderate anisotropy: $`\kappa =50`$, $`\gamma =1/5`$

These values for $`\gamma `$ are typical of YBCO (moderate anisotropy) and BSCCO (strong anisotropy). The parameter $`\kappa `$ for moderate anisotropy is typical of YBCO but for strong anisotropy is about five times smaller than those typical of BSCCO . However, this difference does not alter significantly the phase diagram because only the self-energy depends on $`\kappa `$.

The main results of Ref., summarized in Fig. 3 are as follows. For both moderate and strong anisotropy and for $`\alpha >\alpha _c`$ the vortex lines are nearly parallel to the a-b plane ($`\theta 90^o`$). For strong anisotropy $`\alpha _c=7.9^o`$ and for moderate anisotropy $`\alpha 20^o`$. For strong anisotropy and $`\alpha <\alpha _c`$ the vortex phase at $`H=H_{c1}(\alpha )`$ is a DVLC. For $`\alpha >\alpha _c`$, on the other hand, no significant difference is found between the DVLC and DVL phases, because the vortex lines are tilted at $`\theta 90^o`$, and at such tilt angles the equilibrium chain period is very large. Phase coexistence is also found for strong anisotropy at $`\alpha =\alpha _c=7.9^o`$, but the coexisting phases are a DVLC and DVL, instead of the two DVL found in Ref.. A discontinuous jump in $`\theta `$ takes place as $`\alpha `$ crosses $`\alpha _c=7.9^o`$. For moderate anisotropy the DVL and DVLC phases give nearly identical values for $`H_{c1}(\alpha )`$, so that no conclusion on which is the equilibrium state is reached.

ii) $`\alpha =0^o`$ $`and`$ $`\alpha =90^o`$. In both cases $`\theta =\alpha `$ and the lower critical field is related to the self energy by the usual formula $`H_{c1}(\alpha =0^o,90^o)=4\pi \epsilon _{\mathrm{sf}}(\theta =0^o,90^o)/\mathrm{\Phi }_0`$. For $`\alpha =\theta =0^o`$ the equilibrium VLL is the familiar isotropic superconductor triangular lattice, with $`L_1=L_2=L_\mathrm{\Delta }(\mathrm{\Phi }_0/B\mathrm{sin}60^o)^{1/2}`$, $`\phi =60^o`$ and $`\psi `$ undetermined. For $`\alpha =\theta =90^o`$ the x-direction coincides with the negative $`c`$-axis. The equilibrium VLL is related to that for $`\alpha =\theta =0^o`$ (assuming $`𝐋_1`$ along the x-direction or $`\psi =0^o`$) by the following scaling relations

$`𝐋_1`$ $`=`$ $`L_0\gamma ^{1/2}\widehat{𝐱}`$ (9)
$`𝐋_2`$ $`=`$ $`L_0[\gamma ^{1/2}\mathrm{cos}\phi \widehat{𝐱}+\gamma ^{1/2}\mathrm{sin}\phi \widehat{𝐲}],`$ (10)
where $`L_0=(\mathrm{\Phi }_0/B\mathrm{sin}\phi )^{1/2}`$, with $`\phi =60^o`$ for the triangular lattice. These results also apply if the $`\alpha =0^o`$ VLL is a square lattice, which will be considered in Sec. III, in which case $`\phi =90^o`$ in Eq. (10). The unit cell defined by Eqs. (10) results from compressing the $`\alpha =\theta =0^o`$ VLL unit cell (with $`𝐋_1`$ along the negative $`c`$-direction) by $`\gamma ^{1/2}`$ along the $`c`$-direction and stretching it by $`\gamma ^{1/2}`$ in the ab plane. This transformation conserves the unit cell area $`A_c=L_0^2\mathrm{sin}\phi `$. The $`B`$ vs. $`H`$ curves for $`\theta =\alpha =0^o`$ and $`\theta =\alpha =90^o`$ are also related by scaling as follows

$$B_{90^o}(\frac{HH_{c1}(90^o)}{H_{c1}(0^o)})=\gamma B_{0^o}(\frac{HH_{c1}(0^o)}{\gamma H_{c1}(0^o)}).$$
(11)
iii) High fields ($`HH_{c1}(\alpha )`$). To a good degree of approximation the equilibrium VLL is related to the $`\alpha =\theta =0^o`$ VLL by the scaling relations Ref.

$`𝐋_1`$ $`=`$ $`L_0({\displaystyle \frac{\lambda _\theta }{\lambda _c}})^{1/2}\widehat{𝐱}`$ (12)
$`𝐋_2`$ $`=`$ $`L_0[({\displaystyle \frac{\lambda _\theta }{\lambda _c}})^{1/2}\mathrm{cos}\phi \widehat{𝐱}+({\displaystyle \frac{\lambda _c}{\lambda _\theta }})^{1/2}\mathrm{sin}\phi \widehat{𝐲}].`$ (13)
In this case $`𝐁𝐇`$, with $`B=\mathrm{\Phi }_0/L_\mathrm{\Delta }^2\mathrm{sin}\phi `$, with $`\phi =60^o`$ for the triangular VLL and $`\phi =90^o`$ for the square VLL.

Apart from these cases, little else is known about the vortex phase diagram. Results for $`H=H_{c1}(\alpha )`$, such as the occurrence of a DVLC as the equilibrium state, phase coexistence, and large differences between $`\theta `$ and $`\alpha `$, suggest a low field ($`HH_{c1}(\alpha )`$) phase diagram that is non-trivial and that differs considerably from the high-field one. It is of interest then to obtain a full picture of the phase diagram. This is our main motivation here.

According to the discussion above, for a given $`𝐇`$, $`G`$ depends on five independent variables, $`\theta ,𝐋_\mathrm{𝟏},𝐋_\mathrm{𝟐}`$, or, equivalently, on $`\theta ,L_1,L_2,\phi ,\psi `$ (Fig. 1.b). It can be shown that the dependence of $`G`$ on the angle of rotation of the VLL with respect to $`𝐁`$, $`\psi `$, is such that $`G`$ has a minimum for $`\psi =0`$. We assume that this is the absolute minimum and, from here on, take $`\psi =0`$ in our calculations. In order to simplify the calculations, we fix the value of the angle $`\phi `$ at $`\phi =90^o`$. Thus, in the results reported here, $`𝐋_\mathrm{𝟏}`$ is along the x-axis and coplanar with $`𝐁`$, $`𝐇`$ and the $`c`$-axis, $`𝐋_\mathrm{𝟐}`$ is along the y-axis, and the VLL unit cell is rectangular. We find that setting $`\phi =90^o`$ does not change the general picture for the phase diagram. With these restrictions $`G`$ depends only on three variables: $`\theta ,L_1,L_2`$. In the remainder of the paper we use instead $`\theta ,L_1,B=\mathrm{\Phi }_0/(L_1L_2)`$ as independent variables.

One novel aspect of our calculation is to account for the effects of interactions between vortex lines in the zero-temperature phase diagram. We find that the simple model described above has a non-trivial low-field phase diagram because attractive interactions between vortex lines give an important contribution to $`\epsilon _{\mathrm{int}}`$ at low fields. The reason is that, as discussed above, the equilibrium VLL unit cell has $`𝐋_\mathrm{𝟏}`$ along the x-direction ($`\psi =0`$). This means that the equilibrium phase consists of periodic vortex-line chains aligned along x, with period $`L_1`$, separated from each other by $`L_2`$. According to , the interaction between a pair of vortex lines in the same chain is attractive if their separation $`x`$ is greater than $`x_m(\theta )`$, and repulsive if $`x<x_m(\theta )`$. Thus, as long as $`L_1`$ is not too small compared to $`x_m(\theta )`$, the vortex-line chains contribution to $`\epsilon _{\mathrm{int}}`$ is predominantly attractive and essentially negative. Interactions between interchain vortex lines are always repulsive, thus giving a positive contribution to $`\epsilon _{\mathrm{int}}`$.

In the DVLC state at $`H=H_{c1}(\alpha )`$, these chains are infinitely far apart ($`L_2=\mathrm{}`$) and $`L_1`$ minimizes $`\epsilon _{\mathrm{int}}`$ for the equilibrium $`\theta `$. In this paper, we obtain the vortex phases at $`H=H_{c1}(\alpha )`$ by extrapolating our results for $`H>H_{c1}(\alpha )`$. For strong anisotropy, our results are in agreement with those of Ref.. For moderate anisotropy, we find that the vortex phase at $`H=H_{c1}(\alpha )`$ is also a DVLC for $`\alpha <\alpha _c`$, a result not obtained in Ref..

Our results for the zero-temperature phase diagram can be summarized in the generic phase diagram show in Fig. 2, consisting of three distinct phases, whose main characteristics are discussed in Sec. I.

### B Numerics

The main difficulty to carry out the minimization of $`G`$, Eq. (5), is to evaluate the lattice sum in the expression for $`\epsilon _{\mathrm{in}}`$, Eq. (7). This can only be done numerically. However, the direct evaluation of this sum is impractical for our purposes, because it converges very slowly. To circumvent this difficulty we use the results obtained by Doria , who showed that the lattice sum in Eq. (7) can be transformed into a series that converges much faster than the direct sum.

Our numerical method uses a simulated annealing algorithm to locate the minima of $`G`$. For a given $`𝐇`$, specified by $`H`$ and $`\alpha `$, the minimization of $`G`$ is carried out in the space spanned by $`\theta ,L_1`$, and $`B`$. Our algorithm starts from a convenient choice for these variables and attempts to change them to new values that decrease $`G`$. Once a set of such values is found, another one is searched by repeating the procedure, and so on. This converges to a set corresponding to a minimum of $`G`$ after several steps. For each new set of $`\theta ,L_1,B`$ we calculate $`\epsilon _{\mathrm{in}}`$, using Doria’s fast convergent series. This allows us to run the simulated annealing for a very large number of steps, and to obtain the minima of $`G`$ with high accuracy, except when $`H`$ is very close to $`H_{c1}(\alpha )`$. In this case our method converges slowly, because the VLL unit cell becomes very large and $`\epsilon _{\mathrm{in}}`$ very small. A similar numerical method is described in Ref.

## III Results

Here we report in detail the results obtained by the numerical method described in Sec. II B for strong and moderate anisotropies. To explore the ($`H,\alpha `$) phase diagram we fix $`\alpha `$ and obtain the minima of $`G`$ for several values of $`H>H_{c1}(\alpha )`$. By repeating this procedure for several $`\alpha `$ we obtain the phase diagram. We check that these minima are absolute ones by running the above described simulated annealing algorithm starting from distinct initial states, and selecting the minima with the smallest $`G`$.

### A Strong anisotropy

The main results are summarized in Figs. 4 and 5.

First, we argue that our results are consistent with the prediction, discussed in Sec. I, that at $`H=H_{c1}(\alpha )`$ the vortex configuration for $`\alpha <\alpha _c=7.9^o`$ is a DVLC.

The existence of a DVLC at $`H=H_{c1}(\alpha )`$ for $`\alpha <\alpha _c`$ means that as $`H/H_{c1}(\alpha )1`$ the primitive unit cell sides behave as $`L_1L_{ch}`$(= chain period), $`L_2\mathrm{}`$, and that $`\epsilon _{in}`$, Eq. (5), approaches a negative value, equal to the interaction energy per vortex line in an isolated chain. Our results for $`\alpha =7^o`$ show exactly this behavior. The $`L_1/L_2`$ curve in Fig. 4.c decreases continuously as $`H/H_{c1}(\alpha )1`$, indicating that $`L_1/L_2`$ extrapolates to $`L_1/L_2=0`$ at $`H=H_{c1}(\alpha )`$. On the other hand, the $`L_1`$ curve in (Fig. 5.b) approaches a finite value . We find that the extrapolated values of $`L_1`$, $`\theta `$ and $`\epsilon _{in}`$ at $`H=H_{c1}(\alpha )`$ are $`L_1/\lambda _{ab}=2.8`$, $`\theta =50^o`$ and $`\epsilon _{in}/\epsilon _0=0.04`$. These values agree, within numerical errors, with the chain period, tilt angle and interaction energy per vortex line of the isolated vortex-line chain obtained in Ref..

Next we consider the behavior for $`\alpha <\alpha _c`$ and $`H>H_{c1}(\alpha )`$.

For $`\alpha =7^o`$ our results for the $`L_1`$ vs. $`H`$ curve (Fig. 5.b) show that there are two distinct behaviors. One for $`1<H/H_{c1}(7^o)<1.03`$, where $`L_1`$ increases with $`H`$, and another for $`H/H_{c1}(7^o)>1.03`$ where $`L_1`$ decreases with $`H`$. The $`\epsilon _{in}`$ vs. $`H`$ curve changes slope at $`H/H_{c1}(7^o)=1.03`$, and shows that $`\epsilon _{in}`$ increases continuously with $`H`$, being negative for $`H/H_{c1}(7^o)<1.02`$. A decrease of $`L_1`$ with $`H`$ is what is expected if interactions between vortex lines are purely repulsive. This behavior is identified with phase-III (Fig. 2) and corresponds to the region where $`L_1/L_2>0.6`$ , $`\theta <23^o`$ and $`\epsilon _{in}>0`$. The increase of $`L_1`$ with $`H`$, on the other hand, arises from attractive interactions between vortex lines. This can be seen by plotting $`L_1`$ vs. $`\theta `$ for $`\alpha =7^o`$, as shown in Fig. 5.c. This figure shows that in this range of $`H`$ values, $`L_1`$ is nearly equal to the period of an isolated vortex-line chain with the same $`\theta `$ and has essentially the same $`H`$ dependence. This is only possible if the attractive intrachain vortex-line interactions are strong compared to interchain repulsions. Further evidence for this is found in the negative values for $`\epsilon _{in}`$. We identify the region $`1<H/H_{c1}(7^o)<1.03`$ as phase-I and that for $`H/H_{c1}(7^o)>1.03`$ as phase-III. For $`\alpha =7^o`$ we also find that $`L_2>L_1`$, and that $`L_2`$ decreases with $`H`$ in such a way that $`L_1/L_2`$ increases with $`H`$ as shown in Fig. 4.c.

A similar behavior is expected for other values of $`\alpha <\alpha _c`$. For $`\alpha =5^o`$ we find that the $`L_1/L_2`$ curve shows a downturn as $`H/H_{c1}(\alpha =5^o)`$ decreases (Fig. 4.c), similar to that for $`\alpha =7^o`$, suggesting that the limit value as $`H/H_{c1}(\alpha )1`$ is small. We also find that the limit value of $`\theta `$ as $`H/H_{c1}(\alpha )1`$ (Fig. 4.a) agrees well with that obtained in Ref.. Although, for $`\alpha 5^o`$, we are unable to observe a region where $`L_1`$ increases with $`H`$ and $`\epsilon _{in}<0`$, we believe that the phase diagram for $`\alpha =5^o`$ and $`\alpha =7^o`$ are similar. One difficulty with values of $`\alpha `$ close to $`0^o`$ is that the region where phase-I exists is so close to $`H/H_{c1}(\alpha )=1`$ that the accuracy of our numerical method is insufficient to observe it.

Next we discuss the behavior for $`\alpha >\alpha _c`$ and $`H>H_{c1}(\alpha )`$.

Our results show that $`\theta `$ remains close to $`90^o`$ in the range $`H_{c1}(\alpha )<H<H^{}(\alpha )`$ (Fig. 4.a). We find that for the results in this range the scaling relations discussed in Sec. I. For strong anisotropy $`\gamma =1/\sqrt{200}`$, so that $`\mathrm{tan}\theta >1/\gamma `$ for $`\theta >86^o`$. Accordingly, the field $`H^{}(\alpha )`$ is obtained from the $`\theta `$ vs. $`H`$ curve as the largest $`H`$ for which $`\theta >86^o`$, as indicated in Fig. 4.a. For $`\alpha 10^o`$ and $`H_{c1}(\alpha )<H<H^{}(\alpha )`$ the $`L_1/L_2`$ and $`B`$ vs. $`H`$ curves in Fig. 4 essentially coincide. The common curves are $`L_1/L_2=\gamma =0.07`$ and the $`B`$ vs. $`H`$ curve for $`\alpha =\theta =90^o`$, obtained from that for $`\alpha =\theta =0^o`$ by the scaling relations Eq. (11). The behavior of $`L_1`$, $`L_2`$ and $`\epsilon _{in}`$ in this region is that expected for repulsive interactions: as $`H`$ increases, $`L_1`$ and $`L_2`$ decrease, and $`\epsilon _{in}`$ increase (Fig. 4). From these results we identify the region $`H_{c1}(\alpha )<H<H^{}(\alpha )`$ as phase-II.

For $`H>H^{}(\alpha )`$ and for $`\alpha =10^o`$, $`15^o`$ and $`20^o`$ the tilt angle $`\theta `$ decreases smoothly with increasing $`H`$ (Fig. 4.a). The behaviors of $`L_1`$ and $`\epsilon _{in}`$ with $`H`$ are not smooth. As shown in Fig. 5 there are regions where $`L_1`$ increases with $`H`$ and the $`\epsilon _{in}`$ curve shows a dip, reaching negative values for $`\alpha =10^o`$. This behavior of $`L_1`$ is similar to that for $`\alpha =7^o<\alpha _c`$ in phase-I region. We also find that for $`\alpha 20^o`$ the VLL period along x, $`L_1`$, is close that for an isolated vortex-line chain with the same $`\theta `$, and has a similar $`H`$ dependence (Fig. 5.c). In view of these similarities we identify the region where $`L_1`$ increase with $`H`$ as phase-I. Thus the portion of the phase diagram occupied by phase-I extends into the region $`\alpha >\alpha _c`$ for $`H>H^{}(\alpha )`$, as shown in Fig. 2.

Our results indicate that for strong anisotropy within phase-I the vortex-line chains along x have essentially the same period as an isolated chain tilted at the equilibrium $`\theta `$.

We find that as $`\alpha `$ increases above $`\alpha =20^o`$ the region where $`L_1`$ increases with $`H`$ becomes smaller and the dip in $`\epsilon _{in}`$ decreases. Eventually these features disappear altogether at some $`\alpha `$, so that the phase-I region in the phase diagram is bounded as shown in Fig. 2.

Our results indicate that the transition from phase-II to phase-I is a smooth crossover, with an intermediate region starting just above $`H^{}(\alpha )`$ and ending when $`\theta 70^o`$. In this region the $`L_1/L_2`$ and $`B`$ curves for $`\alpha =10^o`$, $`15^o`$ and $`20^o`$ show only small deviations from the scaling behavior of phase-II, with $`L_1/L_2`$ smaller than the scaling value and the $`B`$ curve slightly above the scaling one (Fig. 4). Behavior characteristic of phase-I, as described above, shows up only for $`\theta 70^o`$. The reason is is that for $`\theta 70^o`$ interactions between intrachain vortex lines are predominantly repulsive because $`L_1x_m(\theta )`$.

For $`\alpha =10^o`$, $`15^o`$ and $`20^o`$ phase-I behavior ends at the $`H`$ value where $`L_1`$ stops increasing with $`H`$. For larger $`H`$, we find that $`L_1`$ decreases smoothly with $`H`$, with $`L_1/L_21`$, and that $`\epsilon _{in}`$ increases with $`H`$ (Fig. 4). We identify these regions as phase-III (Fig. 2).

### B Moderate anisotropy

The results are similar to those for strong anisotropy and are summarized in Figs. 6 and 7.

First we argue that as $`H/H_{c1}(\alpha )1`$ the equilibrium state approaches a DVLC for $`\alpha 20^o`$. As seen in Fig. 6.c, for $`\alpha =10^o`$, $`15^o`$ and $`20^o`$, $`L_1/L_2`$ approaches a small value as $`H/H_{c1}(\alpha )1`$, which we interpret as indicating that $`L_2\mathrm{}`$. The $`L_1`$ vs $`H`$ and $`L_1`$ vs. $`\theta `$ curves, shown in Figs. 7.b and 7.c, suggest that, as $`H/H_{c1}(\alpha )1`$, $`L_1`$ approaches a finite value that is close to the period of an isolated vortex-line chain with the same $`\theta `$ . The $`\epsilon _{in}`$ vs. $`H`$ curve shown in Fig. 7.a becomes negative as $`H/H_{c1}(\alpha )1`$ for $`\alpha =10^o`$ and $`13^o`$ and extrapolates to a negative value at $`H=H_{c1}(\alpha )`$ for $`\alpha =15^o`$ and $`20^o`$. This behavior is only possible for a DVLC state at $`H=H_{c1}(\alpha )`$.

Next we consider the behavior for $`\alpha <\alpha _c`$ and $`H>H_{c1}(\alpha )`$.

For $`\alpha =10^o`$, $`13^o`$ and $`15^o`$ the $`L_1`$ vs. $`H`$ curves in Fig. 7.b show non-monotonic behavior. As $`H`$ increases $`L_1`$ first decreases, then increases and, finally, decreases monotonically. This behavior suggests that attractive interactions between intrachain vortex lines are competing with repulsive interchain ones. This is further supported by the negative $`\epsilon _{in}`$, for the same $`\alpha `$ values, shown in Fig. 7.b. We interpret as phase-I the $`H`$ range starting from $`H=H_{c1}(\alpha )`$ and ending at the $`H`$ value in Fig. 7.b where $`L_1`$ stops increasing. However, the behavior in phase-I for moderate anisotropy differs from that for strong anisotropy. The $`L_1`$ vs. $`\theta `$ curve in Fig. 7.c shows that $`L_1`$ differs considerably from the equilibrium period of an isolated chain with the same $`\theta `$ and has a different $`H`$ dependence.

We find that for $`\alpha >\alpha _c`$ our results agree with phase-II behavior for $`H_{c1}(\alpha )<H<H^{}(\alpha )`$. For moderate anisotropy $`\gamma =1/5`$, and $`\mathrm{tan}\theta >1/\gamma `$ for $`\theta >80^o`$. The field $`H^{}(\alpha )`$ is obtained from the $`\theta `$ vs. $`H`$ curve as the largest $`H`$ for which $`\theta >80^o`$, as indicated in Fig. 6.a. We find that for $`H_{c1}(\alpha )<H<H^{}(\alpha )`$ the $`L_1/L_2`$ and $`B`$ vs. $`H`$ curves shown in Fig. 6 agree well with the scaling predictions Sec. I.

Thus we conclude that for moderate anisotropy the phase diagram is also like that shown in Fig. 2.

## IV Conclusions

In conclusion, we obtain the zero-temperature phase diagram for superconductors with anisotropy parameters typical of BSCCO and YBCO. Three distict phases are found and characterized in detail. Phase-I is where we find new results. The most significant ones are: i)in the case of strong anisotropy we determine how the vortex-line-chain state that occurs at $`H=H_{c1}(\alpha )`$ ($`\alpha \alpha _c`$) is modified by interchain interactions. ii) In the case of moderate anisotropy we show that the vortex state at $`H=H_{c1}(\alpha )`$ is a vortex-line chain state.

In phase-II, we recover the results of Buzdin and Simonov that reveal a large region of the phase diagram over which the vortex lines remain nearly parallel to the ab-plane. Moreover we find the curve, $`H=H^{}(\alpha )`$, around which the vortex lines tilt away from the ab-plane. We also find that this line starts at $`\alpha \alpha _c`$ and $`H=H_{c1}(\alpha _c)`$ and that, as this line is being crossed, the VLL undergoes a structural change that becomes more and more abrupt as $`\alpha \alpha _c`$.

Our results are obtained with the simplifying assumption that $`\phi =90^o`$. We have also carried out the minimization of $`G`$, Eq. (5), with $`\phi `$ as independent variable. The results are identical to those reported above, except for the VLL unit cell shape which, in general, is not rectangular.

Our present results seem to indicate that the equilibrium VLL changes smoothly with $`H`$ and $`\alpha `$. In this case the phase boundaries in Fig. 2 are crossover lines. However, our data for strong anisotropy shows abrupt changes in the slope of the curves for $`\theta `$, $`B`$ and $`L_1`$ vs. $`H`$ for $`H`$ close to $`H_{c1}(\alpha )`$. These could result from a rapid crossover or from true discontinuities. It is found in Ref. that, for strong anisotropy, if $`\alpha `$ is changed and $`H`$ is kept along the $`H/H_{c1}(\alpha )=1`$ curve, there is a discontinuous change in $`\theta `$ at $`\alpha =\alpha _c=7.9^o`$, associated with coexistence of a DVLC state with a DVL one. This coexistence might remain for $`H/H_{c1}(\alpha )>1`$ in the vicinity of $`\alpha _c`$. In the present study we do not attempt to verify if there are true discontinuities. However, our numerical method is capable of doing so. Work on this topic is under way and will be reported elsewhere.

###### Acknowledgements.

This work was supported in part by MCT/CNPq, FAPERJ, CAPES and DAAD. The authors thank Dr. Ernst Helmut Bandt for helpful discussions.

##

In this appendix the formulas that we use to calculate $`\epsilon _{in}(\theta ,L_1,L_2,\phi )`$ numerically are summarized.

Doria showed that the lattice sum in Eq. (7) (with $`𝐋_\mathrm{𝟏}`$ along x) can be transformed in a fast convergent series . The results are:

$`\epsilon _{in}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Phi }_0}{4\pi \lambda _{ab}}}\right)^2\{{\displaystyle \frac{\lambda _\theta }{\lambda _c}}V_0(0)(1\gamma ^2)\mathrm{cos}^2\theta {\displaystyle _0^1}{\displaystyle \frac{du}{\sqrt{c_1(u)c_2(u)}}}[\mu _0^2(u){\displaystyle \frac{V_0(u)}{\mu _0^2(u)}}]+`$ (14)
$`+`$ $`{\displaystyle \frac{|\mathrm{cos}\theta |}{\gamma ^2}}\mathrm{ln}[{\displaystyle \frac{\lambda _\theta +\lambda _c|\mathrm{cos}\theta |}{\lambda _{ab}(1+|\mathrm{cos}\theta |)}}]\}`$ (15)
where $`V_0`$ and its derivative are given by,

$`V_0(u)`$ $`=`$ $`{\displaystyle \frac{1}{2\mu _0\mathrm{tanh}\left(\mu _0\sigma /2\right)}}+2{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{s=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{cos}(m\chi s)}{\sqrt{\mu _0^2+m^2}}}\mathrm{exp}(\sqrt{\mu _0^2+m^2}\sigma s)+`$ (16)
$`+`$ $`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{\sqrt{\mu _0^2+m^2}}}{\displaystyle \frac{1}{m}}\right)+\mathrm{ln}(\mu _0/2)+Ec`$ (17)
$`\mu _0^2(u){\displaystyle \frac{V_0(u)}{\mu _0^2(u)}}`$ $`=`$ $`{\displaystyle \frac{1}{4\mu _0}}[{\displaystyle \frac{1}{\mathrm{tanh}\left(\mu _0\sigma /2\right)}}+{\displaystyle \frac{\mu _0\sigma /2}{\mathrm{sinh}^2(\mu _0\sigma /2)}}]+`$ (19)
$`+`$ $`\mu _0^2{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{s=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mathrm{cos}(m\chi s)}{\mu _0^2+m^2}}\mathrm{exp}(\sqrt{\mu _0^2+m^2}\sigma s)[{\displaystyle \frac{1}{\sqrt{\mu _0^2+m^2}}}+\sigma s]+`$ (20)
$`+`$ $`{\displaystyle \frac{\mu _0^2}{2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}({\displaystyle \frac{1}{(\mu _0^2+m^2)^{3/2}}}`$ (21)
where $`Ec=0.5772\mathrm{}`$ is Euler’s constant and $`\chi `$ is defined as $`\chi =(2\pi L_2/L_1)\mathrm{cos}\phi `$. The functions of $`c_1`$, $`c_2`$, $`\sigma `$ and $`\mu _0`$ are defined as:

$`c_1(u)={\displaystyle \frac{\lambda _\theta ^2}{\lambda _c^2}}u(1\gamma ^2)\mathrm{cos}^2\theta ,`$

$`c_2(u)=1u(1\gamma ^2),`$

$`\sigma =2\pi {\displaystyle \frac{L_2}{L_1}}\sqrt{{\displaystyle \frac{c_1(u)}{c_2(u)}}}\mathrm{sin}\phi ,`$and

$`\mu _0={\displaystyle \frac{L_1}{2\pi \lambda _c\sqrt{c_1(u)}}}.`$

In the expression for $`\epsilon _{in}`$, Eq. (15), $`V_0(0)`$ stands for $`V_0(u)`$, Eq. (17), for $`u=0`$, or, in other words with the functions $`c_1`$, $`c_2`$, $`\sigma `$ and $`\mu _0`$ in this equation evaluated at $`u=0`$.

In our numerical calculations of $`\epsilon _{in}`$ we evaluate $`V_0(0)`$ by using the following generic approximation for the summations in eqs.(17) and (21)

$$\underset{n=1}{\overset{\mathrm{}}{}}g(n)\underset{n=1}{\overset{M}{}}g(n)+\frac{1}{2}_M^{\mathrm{}}𝑑x(g(x)+g(x+1)),$$
(22)
where $`g(x)`$ is an analytic function. By using $`M=200`$ we obtain excellent numerical precision with modest computer run times.

The integral in eq.(15) is calculated using standard integration algorithms found in the literature, up to a precision of one part in $`10^8`$. We run typical calculations of 10-20,000 Metropolis steps for each run. The results obtained are easily reproducible.


text:
# Characteristic lengths at moving contact lines for a perfectly wetting fluid: the influence of speed on the dynamic contact angle

## 1 Introduction

One area of fluid mechanics that has been the subject of a large admixture of analysis, experiment and speculation is the subject of the moving contact line. A typical situation, common in many coating processes, refers to the contact line at the intersection of solid, liquid and gas regions, where the three-phase line moves relative to a solid substrate. A basic research question in this subject stems from the violation of the no-slip condition in the immediate neighborhood of the three-phase line of contact (e.g. Huh & Scriven 1971, for reviews see Dussan V. (1979), de Gennes (1985), and Kistler (1993)). As a result, within the usual continuum analysis, the stress diverges as the contact line is approached and the energy per unit length of the moving contact line is unbounded. This result may be viewed as an embarrassment of continuum modeling, but, in fact, it does indicate the need for a small cut-off length scale in macroscopic theories, as well as some more input from the physics at smaller length scales to properly interpret the meaning of any such cut-off scale.

Perhaps the most basic feature of this problem is the aim to relate the local dynamic contact angle $`\theta _d(x)`$, which is the arc tangent of the slope of the interface at a distance $`x`$ from the contact line, to the local speed $`U`$ with which the contact line moves over the substrate. For the case of a perfectly wetting fluid (vanishing equilibrium contact angle $`\theta _{eq}=0`$) and small $`\theta _d`$, one finds $`\theta _d^3U`$, which is known as Tanner’s law. Theoretical justification for this result has been given (e.g. de Gennes 1985) and various generalizations have been offered. In dimensionless form, the speed is reported in terms of the capillary number $`𝒞=U\eta /\gamma `$, which measures the relative importance of viscous to surface tension forces, where $`\eta `$ is the fluid viscosity and $`\gamma `$ the interfacial tension. In fact, the functional form for the contact angle-speed relation is commonly written for small angles as $`\theta _d^3(x)9𝒞\mathrm{ln}(x/\mathrm{}_{micro})`$, where $`\mathrm{}_{micro}`$ is generally taken as a molecular length (e.g. Leger & Joanny 1992). Typically the capillary number varies over many orders of magnitude; values $`10^7<𝒞<10^1`$ are common. The prefactor in this formula can be important for interpreting experimental data, and so it is reasonable to interrogate more closely the functional dependence on speed.

In this communication we wish to comment on one aspect of the moving contact line problem that has, perhaps surprisingly, been largely neglected and/or unappreciated. In particular, we note that detailed models for the perfectly wetting situation actually yield a dynamic contact angle versus speed relation

$$\theta _d^3(x)9𝒞\mathrm{ln}\left(\frac{x}{\mathrm{}_{micro}}𝒞^\beta \right),$$
(1)
where $`\beta `$ depends on the physical model introduced in the neighborhood of the contact line. We do not believe that it is necessarily appropriate to simply suppress the additional dependence on speed (i.e. $`𝒞`$) by replacing the argument of the logarithm by either $`\mathrm{}_{macro}/\mathrm{}_{micro}`$, where these two lengths scales are taken as constants, or $`x/\mathrm{}_{micro}`$. Because of the large variation in $`𝒞`$, not including this additional factor of capillary number when using (1) to interpret dynamical experiments may lead to significant discrepancies between theory and experiment. Here we outline the basic idea behind (1) and present experimental evidence that supports the above interpretation.

Another point that has received insufficient attention is the range of validity of equation (1). Near the contact line, (1) breaks down where $`x`$ is of the same order as $`\mathrm{}_{micro}`$. This restriction is evident as the general structure comes from a balance of viscous and surface tension forces alone. Not surprisingly, we estimate below that the microscopic scale is between several Angstroms and tens of Angstroms, depending on the microscopic forces assumed to be acting near the contact line.

Towards large scales, $`x`$ is commonly taken to be a static scale such as the capillary length or the size of a spreading drop (de Gennes 1985). Nevertheless, it should be noted that the flow near a moving contact line often resembles a coating flow, similar to the classical problem studied by Landau and Levich (e.g. Levich 1962). This leads to the appearance of another, dynamical length scale, that can become much smaller than the capillary length as the capillary number is small, which is typically the case. Thus a meaningful comparison between (1) and a macroscopic measurement of the dynamical contact angle might require a spatial resolution significantly below 1/10 or even 1/100 of the capillary length.

In the next section we will introduce two different models commonly used to treat moving contact line problems such as a spreading drop or a tape plunging into a pool of fluid. Then, in the third section we show that the lubrication equations corresponding to both models have similarity solutions for the interfacial shape that fix the functional dependence on the capillary number. The fourth section discusses the dynamical problem that equation (1) has to be matched to on an appropriate outer length scale. In the fifth section we explain measurable consequences of the two models for the dynamic contact angle, and discuss an experiment that helps to distinguish between them. We close with a summary and possible directions of future work.

## 2 The model

The usual dynamic balance for the steady flow “far” from the contact line involves capillary and viscous stresses. As the contact line is approached, the capillary-viscous flow leads to a stress singularity. A number of different physical effects have been suggested to relieve the singularity, and these either account for the fact that on very small length scales van der Waals forces act to maintain a finite thickness liquid layer on the solid substrate, or that at very high shear rates the boundary conditions and the transport coefficients of the fluid are likely to be altered. Which model is appropriate might depend on the physical system at hand, or be a combination of the above. In Table 1 we provide a short overview of proposed physical models for flow in the neighborhood of a contact line; see also McKinley & Ovryn (1998). Most recently, there has been a considerable effort to base the understanding of the contact line physics on a microscopic, particle-based description, see for example Koplik et al (1989), Ruijter et al. (1999), and Abraham et al. (2002). The so-called “diffuse interface model” (see e.g. Seppecher (1996), Chen et al. (2000), and Pomeau (2002)) represents an intermediate approach, which models the liquid-gas interface as a Cahn-Hilliard fluid. This allows for example the extraction of effective interface equations (Pismen & Pomeau (2000)) different from those proposed by de Gennes’s (1985).

Below we restrict our attention to two different models which have proved particularly popular. The results are sufficient to highlight the measurable differences between different physical mechanisms. In model I, due to de Gennes and coworkers (e.g. Hervet & de Gennes 1982, de Gennes 1985), van der Waals forces are taken into account, so very close to the contact line there is a balance between surface tension and van der Waals stresses alone. In model II, proposed for example by Huh and Mason (1977) and Hocking (1977), the fluid is allowed to slip across the solid surface over a small slip length.

For simplicity, we only consider the case of perfectly wetting fluids, i.e. of zero equilibrium contact angle. Consistent with the local balances, the interface near the contact line remains nearly flat and we can use lubrication theory to describe the fluid motion. This approach amounts to a significant simplification of the mathematical treatment relative to the full two-dimensional flow problem (Cox 1986), but agrees with the full calculation when the dynamic contact angle is small. There are numerous indications that the small-angle theory in fact remains valid for slopes of order unity. For example, $`\theta _d^3/9`$ in equation (1) differs by only $`2\%`$ from the full expression (Cox 1986), derived without the benefit of lubrication theory, up to a slope of 1. Also, Fermigier & Jenffer (1991) reported that small-angle theory holds experimentally up to an angle of $`100^{}`$.

To be able to describe an experiment like a flat plate plunging with velocity $`U`$ into a reservoir of fluid (see Figure 1), it is necessary to include other terms beyond the lubrication terms close to the contact line. Namely, we keep the full curvature term and include gravity. By doing this, the model is able to describe the crossover to a purely static, horizontal surface far away from the dynamical region. A basic unit of length is the capillary length $`\mathrm{}_c=\sqrt{\frac{\gamma }{g\rho }}`$, which dictates the scale of the interface curvature far away from the contact line. In the van der Waals model I, which accounts for pressure variations owing to capillary, van der Waals, and gravitational forces, the equation for the stationary profile $`h(x)`$ (cf Figure 1) is

$$\frac{3𝒞}{h^2}=\kappa ^{}+3a^2\frac{h^{}}{h^4}\mathrm{}_c^2,$$
(2)
where $`\kappa `$ is the curvature and a prime refers to differentiation with respect to $`x`$. Note that a positive $`𝒞`$ corresponds to the plate plunging into the fluid. A brief derivation of model I, as well as model II below, is given in appendix 1. The microscopic length parameter $`a`$, defined by

$$a^2=\frac{A}{6\pi \gamma },$$
(3)
measures the strength of van der Waals forces relative to interfacial forces and is typically very small (on the order of Angstroms).

Another distinct approach for the flow near the contact line is to introduce slip at the boundary, consistent with allowing the contact line to move parallel to the wall at a finite speed; the slip is a function of the shear rate. The simplest such law, introduced by Navier in the same paper that also enunciated the Navier-Stokes equation (Navier 1823), is

$$u|_{y=0}U=\lambda \frac{u}{y}|_{y=0}$$
(4)
(see also Huh & Scriven 1971). Here $`U`$ is the speed of the moving boundary, $`y=0`$ denotes the solid-liquid boundary, and $`\lambda `$ is the so-called slip length. A more complicated version of (4), in which $`\lambda `$ is itself a non-linear function of the shear rate, has been proposed in Thompson & Troian (1997). A standard calculation (appendix 1), leads to the analogue of equation (2) for the slip model II,

$$\frac{3𝒞}{h^2}=\kappa ^{}3\lambda \frac{(1\kappa ^{})}{h}\mathrm{}_c^2.$$
(5)
The slip length $`\lambda `$ is usually considered to be on the order of tens of Angstroms.

## 3 Scaling solutions near the contact line

We now focus on the immediate neighborhood of the contact line, which we assume to be at $`x=0`$. Owing to the flatness of the interface we can assume that $`\kappa h^{\prime \prime }`$ and gravitational influences can be neglected, but dynamical (viscous) effects have to be included. In the case of model I, equation (2) reduces to

$$\frac{3𝒞}{h^2}=h^{\prime \prime \prime }+3a^2\frac{h^{}}{h^4}.$$
(6)
To make the dependence on parameters explicit, we note that (6) has the exact scaling solution

$$h(x)=\frac{a}{𝒞^{1/3}}\varphi _1(a^1𝒞^{2/3}x),$$
(7)
where $`\varphi _1`$ depends on the similarity variable $`\xi _1=a^1𝒞^{2/3}x`$ and satisfies the equation

$$\frac{3}{\varphi _1^2}=\varphi _1^{\prime \prime \prime }+3\frac{\varphi _1^{}}{\varphi _1^4}.$$
(8)
Similarly, the lubrication approximation for model II gives

$$\frac{3𝒞}{h^2}=h^{\prime \prime \prime }\frac{3\lambda }{h}h^{\prime \prime \prime }.$$
(9)
In this case the scaling solutions are

$$h(x)=\lambda \varphi _2(\lambda ^1𝒞^{1/3}x),$$
(10)
where the similarity variable is now $`\xi _2=\lambda ^1𝒞^{1/3}x`$ and the similarity equation is

$$\frac{3}{\varphi _2}=\varphi _2^{\prime \prime \prime }\frac{3}{\varphi _2}\varphi _2^{\prime \prime \prime }.$$
(11)
Far away from the contact line in units of the microscopic lengths $`a`$ and $`\lambda `$, respectively, the solutions should be the same, resulting from a balance of classical viscous forces and surface tension. Indeed, as $`\xi \mathrm{}`$ one finds to leading order

$$\varphi _{1,2}(\xi )3^{2/3}\xi \left[\mathrm{ln}(\xi b_{1,2})\right]^{1/3},\xi 1,$$
(12)
where the numerical constants $`b_{1,2}`$ have to be determined by numerical integration starting from the contact line.

The boundary condition at the contact line incorporates the wetting behavior of the fluid. The basic assumption is that there is a microscopic length scale on which static forces dominate over dynamical ones, and a static profile can be assumed microscopically close to the contact line.

In model I, following Hervet & de Gennes (1984), we are going to match $`\varphi _1(x)`$ to a “maximal” film solution, corresponding to very strong wetting, whose thickness only goes to zero at (minus) infinity. However this maximal solution very closely approximates parabolic solutions of (6) that go to zero at some finite contact line position (Hervet & de Gennes 1984). To leading order, we take the film solution to be of the form

$$\varphi _1(\xi _1)=\frac{1}{\xi _1}+ϵ\mathrm{exp}\left\{\xi _1^3/\sqrt{3}\right\},$$
(13)
and defer further details to appendix 2. Using (13) as an initial condition with adjustable parameter $`ϵ`$, we integrate (8) towards $`\xi _1\mathrm{}`$. The parameter $`ϵ`$ is fixed to select the solution with vanishing curvature at infinity. Figure 2 compares this solution with the asymptotic form (12). We plot the rescaled slope $`\varphi _1^{}`$ from the solution of (8) as the full curve and equation (12) with $`b_1=1.44`$ as the dashed curve. This numerical value for $`b_1`$ differs significantly from $`b_1=0.43^{1/6}0.48`$ given in de Gennes (1985), accounting for differences in normalization. We believe the difference is simply due to the large values of $`\xi `$ necessary for integration until a true asymptotic value is reached.

Model II, on the other hand, can be extended down to $`\varphi _2(0)=0`$, since the stress singularity was successfully removed. A vanishing equilibrium contact angle $`\theta _{eq}=0`$ can thus be implemented by taking the boundary condition $`\varphi _2(0)=\varphi _2^{}(0)=0`$ at the contact line and integrating (9) towards $`\xi _2\mathrm{}`$. The corresponding value of $`b_2`$ for the case of the Navier slip law was given in Hocking (1992). Thus the two constants, which establish the form of the interface profile, are

$$b_11.44,b_23^{1/3}\mathrm{exp}(0.74/3)1.85.$$
(14)
## 4 Crossover to Landau-Levich-type behavior

We now estimate the range of validity of the solution (12) as one moves farther away from the contact line. These ideas have close analogy to the classical analysis of Landau and Levich of a dynamical lubrication film (Levich 1962). Evidently, the contact line physics plays no role far away from the contact line, so the relevant lubrication equation is

$$\frac{3𝒞}{h^2}=h^{\prime \prime \prime },$$
(15)
which is to be matched to a static profile at large distances. This equation has the general solution

$$h(x)=\mathrm{}_c𝒞^{\alpha _1+1/3}f(x/(\mathrm{}_c𝒞^{\alpha _1})),$$
(16)
which has to be matched to a static meniscus on the capillary scale. This static solution is characterized by a curvature that is approximately constant, thus $`h^{\prime \prime }(x)`$ must be independent of $`𝒞`$. This fact forces $`\alpha _1=1/3`$, so that we have

$$h(x)=\mathrm{}_c𝒞^{2/3}f(x/(\mathrm{}_c𝒞^{1/3})),$$
(17)
which implies that the crossover will occur on a scale $`\mathrm{}_c𝒞^{1/3}`$, on which the logarithmic dependence in equation (12) for the slope begins to fail.

The crossover to scaling of the form of equation (17) is demonstrated in Figure 2, by showing a full solution of equation (2), rescaled according to (7). Results are given for two (small) values of the capillary number differing by a factor of ten. Again, the free parameter in the maximal film solution (13) is used to shoot for the flat interface corresponding to the surface of the fluid-filled container. For small values of $`\xi _1`$, the solution corresponds to the lubrication form given before, while on a scale $`x/\mathrm{}_c𝒞^{1/3}`$ the transition to the Landau-Levich region is observed. In rescaled coordinates $`\xi _1`$ the location of this crossover should thus be proportional to $`𝒞`$ itself, as is clearly seen from Figure 2. We have chosen the smaller of the two values of $`𝒞`$ such that the region over which the asymptotic form (12) of the interface can be applied is zero, to highlight possible problems in comparing asymptotic solutions with experimental data. To interpret the measured dynamical contact angle equation (12) is no longer sufficient, but the full solution of the similarity equation (8) has to be considered.

## 5 Comparison with experiment

It is common practice in the literature (e.g. de Gennes 1985, Cox 1986) to consider the derivative of the profile $`h(x)`$, evaluate it at some macroscopic distance from the contact line $`x=\mathrm{}_{macro}`$, and to interpret the slope of the interface in terms of the so-called “dynamical contact angle”, $`\mathrm{tan}\theta _d(x)=dh/dx`$. This approach is the common one taken in experiments as well. Thus, using the solution (12) in the similarity forms (7) and (10) for models I and II, respectively, and neglecting lower-order terms, we obtain

$$\theta _{dyn}^3(x)=9𝒞\mathrm{ln}(\mathrm{}_{macro}/L_{1,2}),$$
(18)
where $`L_{1,2}`$ are microscopic lengths appropriate for each model. There are two fundamental issues with this approach: First, depending on the experimental system it is not clear what is the best choice for $`\mathrm{}_{macro}`$. Second, what is usually taken as a fixed microscopic length $`L_{1,2}`$ is actually strongly dependent on the capillary number. Namely, the two models give

$$L_1=a𝒞^{2/3}/b_1,\text{and}L_2=\lambda 𝒞^{1/3}/b_2.$$
(19)
as established in (7) and (10), respectively.

In particular, the $`𝒞`$ dependence that appears in the microscopic length is different in the two models. It is also clear that it is impossible to interpret $`L_{1,2}`$ directly in terms of some fixed microscopic length near the contact line, but rather it is a dynamical quantity. To our knowledge, this fact has never been appreciated in either theoretical or experimental work. This observation appears to be significant, since by comparing the $`𝒞`$-dependence it potentially allows one to distinguish between different microscopic models from a macroscopic measurement. Chen and Wada (1989) imaged the profile near the contact line of a spreading droplet and so provided the first experimental confirmation of (18). However, owing to the small range of capillary numbers studied, it is difficult to distinguish between the two lengths $`L_{1,2}`$ defined in (19). Below we will therefore concentrate on another experiment (Marsh, Garoff & Dussan 1993), which allowed $`𝒞`$ to be varied over more than two orders of magnitude.

It is also interesting to note that the logarithmic dependence on capillary number was only obtained in the fully nonlinear treatment outlined above. In the classical studies of the flow in the neighborhood of the dynamic contact angle, for example in Cox (1986), a matched asymptotic analysis is used which further assumes a form involving integer powers of $`𝒞`$ and does not recognize that the scale of the inner region can itself involve $`𝒞`$. However the form of the asymptotics (18), implying a logarithmic dependence on $`𝒞`$ in the full solution, clearly shows that the profile cannot be expanded in integer powers of C. As a result, the classical analyses are not able to identify the sort of dependencies given by (19). Although these dependencies are only logarithmic, as mentioned above, the capillary number often varies through many orders of magnitude in experiments, so the logarithm in (18) can in general not be approximated by a constant, as is most often done (e.g. King 2001).

In a recent experiment, the effect of large variations of $`𝒞`$ on the contact line was investigated very carefully by Marsh et al. (1993), who measured the dynamic contact angles on a cylinder plunging at an angle into a liquid bath. These authors essentially used the form (18) to fit the whole shape of the interface close to the contact line, and included static contributions to account for the effects of surface tension and gravity away from the contact line. (Note that instead of the third power on the left hand side, they actually used a more complicated function $`g(x)`$, but which becomes $`g(x)x^3/9`$ for small arguments. This limit is relevant for the small angle case we are studying here.) This approach leaves out dynamical effects of the kind predicted by Landau and Levich (Levich 1962), which are important on an intermediate scale between the microscopic ones and the capillary length, and should be taken into account in a more refined theory. Marsh et al. (1993) treat the static contact angle $`\theta _{eq}`$ (called $`\theta _{act}`$ by the authors) as a free, and possibly $`𝒞`$-dependent, parameter, to be determined from experiment. The authors conclude that for their system $`\theta _{eq}=0`$, which is the case treated here.

From the fit of (18) to their data, Marsh et al. (1993) extract a length $`L`$, which is found to depend significantly on capillary number, as suggested by (19)! They also report $`L`$ to be independent of the tilt angle $`\alpha `$ within experimental error, which further emphasizes that the response is dominated by local features. In Figure 3 we present a plot of the measured length as function of capillary number, and compare it with the slopes suggested by the van der Waals model I and the Navier slip model II, respectively. Although it is difficult to draw firm conclusions owing to the large scatter in the data, 2/3 seems to be favored. Using the two different fits plotted in Figure 3, we are also able to determine the cut-off lengths $`a`$ and $`\lambda `$, assuming that the corresponding physical mechanism is really relevant for the particular materials involved. We find $`a4`$Åfor the van der Waals model and $`\lambda 30`$Åfor the slip model. Using the value of $`A=10^{20}\text{J}`$ (Russel, Saville & Schowalter 1989) for the Hamaker constant for water and an adjacent solid surface, and $`\gamma =0.07\text{N/m}`$, we find $`a1`$Å, consistent with the above value. However, it is important to keep in mind that there is no reason why a single cutoff mechanism should necessarily dominate in the experiment, which would lead to still other exponents. Additional mechanisms for relieving the contact line singularity are listed in Table 1; which is the dominant mechanism could also depend on capillary number and, in particular, on the type of solid substrate or fluid involved.

## 6 Conclusions

We have considered flow local to a moving contact line using a lubrication approach. Our basic message has been to indicate that the appearance of logarithmic corrections in capillary number to the usual “Tanner’s law”, $`\theta _d^3𝒞`$, are a general feature of the mechanical response. The interpretation of the results is that the “microscopic” length scale that is involved when supplying a small-scale cut-off to relieve the well-known stress singularity in the moving contact line problem is a dynamical (speed dependent) quantity. Experimental data consistent with this interpretation is cited, and microscopic parameters that come from the comparison with different theories agree with physical considerations.

Since important applications of contact line theories apply to angles up to $`180^{}`$, it would be very useful to extend the lubrication theory considered here to a full two-dimensional treatment of the flow in the corner region. We suspect that large parts of the calculation in Cox (1986), which erroneously assumes a regular expansion in $`𝒞`$, could be adapted to a proper similarity description. This means that the dynamical contact angle has the scaling form $`\theta _d=f(𝒞,x/L_{1,2})`$, where $`L_{1,2}`$ is one of the dynamical length scales defined by (19), with corresponding forms for the velocity field as well. Such a theory might be able to explain more recent experiments (Chen, Ramé & Garoff 1995) on moving contact lines performed at higher capillary numbers, yielding dynamical contact angles of up to $`155^{}`$.

Furthermore, the velocity of the contact line relative to the substrate is in general not perpendicular to the contact line, as highlighted in recent experiments of droplets running down an inclined plane (Podgorski, Flesselles & Limat 2001). In this case the flow is truly three-dimensional, and it may no longer be sufficient to simply project the velocity onto the normal to the contact line (Blake & Ruschak 1979). Such a three-dimensional description would be necessary to complete our understanding of corner singularities that form at the back of running drops (Stone et al. 2002), and may apply to a range of other contact line phenomena as well.

## 7 Appendix 1: Derivation of lubrication equations

Here we briefly recall the derivation of the interface, or lubrication, equations for thin viscous films (Levich (1962)). For pressure-driven flow along the surface and absorbing hydrostatic pressure variations into the pressure $`p`$, the velocity parallel to the plate can be represented as a second-order polynomial

$$u=a_0+a_1y+y^2\frac{p^{}}{2\eta },$$
(20)
where $`y`$ is the distance normal to the plate. At the free surface $`y=h(x)`$, shear gradients $`u/y`$ vanish, giving $`a_1=p^{}h/\eta `$. Finally, from the slip condition (4) we have $`a_0=U\lambda p^{}h/\eta `$.

Since the contact line is stationary, the mass flux through the film is zero everywhere, $`_0^hu(y)dy=0`$, and thus

$$0=(U\lambda p^{}h/\eta )hp^{}h^3/3\eta $$
(21)
is the equation for the film profile $`h(x)`$. In the presence of van der Waals forces, the dynamic pressure in the liquid is

$$p=\gamma \kappa \frac{A}{6\pi h^3}g\rho x,$$
(22)
where $`A`$ is Hamaker’s constant. Substituting (22) into (21) and assuming $`\lambda =0`$ gives (2), while $`A=0`$ at finite $`\lambda `$ leads to (5).

## 8 Appendix 2: expansion for the maximal film

Here we give some more details on the solution of (6) for the “maximal film” of Hervet and de Gennes (1984). The general form of the film profile is

$$\varphi _0=\frac{1}{\xi }\underset{i=0}{\overset{\mathrm{}}{}}\frac{a_i}{\xi ^{6i}},$$
(23)
where we denote the similarity variable by $`\xi `$. This expansion has no free parameters, as the values of the coefficients $`a_i`$ are obtained directly from substituting (23) into (6). We find

$$a_0=1,a_1=2/5,a_2=1764/275,\mathrm{}.$$
(24)
However, there is a one-parameter family of solutions of (6) that decay for $`\xi \mathrm{}`$. This solution is found by linearizing around the base solution (23), i.e. $`\varphi (\xi )=\varphi _0(\xi )+\delta (\xi )`$:

$$\delta (6\varphi _0+4\varphi _0^3\varphi _0^{\prime \prime \prime })3\delta ^{}+\delta ^{\prime \prime \prime }\varphi _0^4=0.$$
(25)
Equation (25) is solved using a WKB-type ansatz,

$$\delta (\xi )=ϵ\mathrm{exp}\left\{\frac{\xi ^3}{\sqrt{3}}+\mathrm{}\right\}.$$
(26)
The $`O(\xi ^0)`$ contribution in the exponent turns out to be a logarithm, so the full structure is

$$\delta =\frac{ϵ}{\xi ^2}\mathrm{exp}\left\{\underset{i=0}{\overset{\mathrm{}}{}}\frac{b_i\xi ^{33i}}{\sqrt{3}}\right\},$$
(27)
and the coefficients are found to be

$$b_0=1,b_2=0,b_3=32/15,b_4=9\sqrt{3}/5,\mathrm{}.$$
(28)
Thus the general form of the solution in the film region is

$$\varphi (\xi )=\varphi _0(\xi )+\delta ,$$
(29)
with a single free parameter $`ϵ`$. An alternative description would be an expansion of the form

$$\varphi (\xi )=\underset{i=1}{\overset{\mathrm{}}{}}\frac{c_i}{\xi ^i},$$
(30)
with $`c_1=1`$ and $`c_2`$ a free parameter. However, the convergence of the asymptotic series (30) turns out to be very bad, as perhaps is to be expected from the structure of the WKB solution.

###### Acknowledgements.

We thank Cyprien Gay, Pirouz Kavehpour, Laurent Limat, Gareth McKinley, Thomas Podgorski and David Quéré for helpful conversations. HAS thanks the Harvard MRSEC for partial support of this research.