text:
# Singular current response from isolated impurities in 𝑑-wave superconductors

## Abstract

The current response of a $`d`$-wave superconductor containing a single impurity is calculated and shown to be singular in the low-temperature limit, leading in the case of strong scattering to a $`1/T`$ term in the penetration depth $`\lambda (T)`$ similar to that induced by Andreev surface bound states. For a small number of such impurities, we argue this low-$`T`$ upturn could be observable in cuprate superconductors.

Introduction. The quasiparticle excitations near the line nodes in the $`d`$-wave superconducting order parameter of the cuprate high-$`T_c`$ materials lead to well-known singularities in the current response of such systems, resulting in nontrivial magnetic field dependences of the London penetration depth in the Meissner stateyipsauls , and of the specific heat in the Abrikosov vortex statevolovik . These excitations are also responsible for the marginal thermodynamic stability of the $`d`$-wave state itself, which has been relatedschopohl ; lietalcomment  to the famous linear-$`T`$ temperature dependence of the penetration depth $`\lambda (T)`$hardy .

Recently, it has been pointed out by Walter et al.walteretal  and Barash et al.barashetal , that, in addition to effects arising from extended quasiparticle states alone, there are $`1/T`$ contributions to the penetration depth below a crossover temperature $`T_{m0}\sqrt{\xi _0/\lambda _0}T_c`$, where $`\xi _0`$ is the coherence length and $`T_c`$ the critical temperature, due to zero-energy surface bound states. These states arise within a semiclassical picture when a quasiparticle reflected from the surface experiences a change of sign in the order parameter; the upturns therefore only occur when the surface normal of the sample makes an angle close to $`\pi /4`$ with the crystal axes and hence the antinodal directions of the $`d`$-wave order parameter. A $`\lambda (T)1/T`$ term in the penetration depth has in fact been reported by Carrington et al.prozorov , and attributed to the paramagnetic current carried by such states which dominates at low $`T`$. The field and interface angle dependence also appear to be in rough agreement with the predictions of Barash et al.barashetal . Upturns correlated with disorder introduced by Zn atoms observed by Bonn et al.bonnetal  seem to vary from sample to sample, and are not currently understood.

Most works in this area have treated disorder within effective-medium approximations which predict broadening of low-energy quasiparticle states by a residual rate $`\gamma `$, which is also roughly proportional to the residual density of states at zero energy $`N(0)`$. Impurities are thus assumed to “smear out” the nodes of the gap, and therefore inevitably cut off the singular behavior. For example, in the work of Barash et al.barashetal , no upturn in $`\lambda (T)`$ at low $`T`$ is observed if $`\gamma >T_{m0}`$. In the impurity-dominated regime, a quadratic temperature dependence is generally to be expectedgrossetal . Impurity physics is therefore thought to competelietalcomment  with nonlocalkosztin  and nonlinearyipsauls  effects in the Meissner state, and $`\lambda (T)`$ changes from linear to quadratic below whichever scale is largest.

There may be reasons, however, to doubt results obtained for the penetration depth using the self-consistent $`T`$-matrix approximation (SCTMA), which defines a translationally invariant effective medium for the impurity-averaged dirty $`d`$-wave superconductor. First, it has been shown to break down in two dimensionstsvelik  and corrections due to weak localization and correlated order parameter response have been under intense investigation recently.pjhjltp  Secondly, as with any effective medium theory, it must break down when impurities are sufficiently isolated. Finally, there are clear indications that single impurity bound states in unconventional superconductors are analogous in some respects to surface Andreev states.goldbart  Here we point out that under some circumstances isolated strong impurities can themselves make singular contributions to the penetration depth and other thermodynamic quantities. While we can present only estimates for the effect of a thermodynamically finite density of dilute impurities, we perform an exact calculation for the current response of a $`d`$-wave superconductor in the presence of a single strong scatterer.

Current response for single impurity. We consider a pure $`d`$ wave superconductor described by Nambu propagator $`\underset{¯}{G}_𝐤^0(\omega )=(i\omega _n\tau _0+\xi _𝐤\tau _3+\mathrm{\Delta }_𝐤\tau _1)/D_𝐤`$, where $`D_𝐤\omega _n^2+\xi _𝐤^2+\mathrm{\Delta }_𝐤^2`$, $`\xi _𝐤ϵ_𝐤\mu `$ is the 1-electron spectrum measured relative to the Fermi level, and the $`\tau _i`$ are Pauli matrices. The $`d`$-wave order parameter on the model cylindrical Fermi surface parametrized by angle $`\varphi `$ is $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0\mathrm{cos}2\varphi `$. In the system with one single $`\delta `$-function impurity of strength $`V_0`$ located at position $`𝐑`$, the Green’s function is exactly given by

$`\underset{¯}{G}_{\mathrm{𝐤𝐤}^{}}(\omega )=\underset{¯}{G}_𝐤^0(\omega )\delta _{\mathrm{𝐤𝐤}^{}}+e^{i(𝐤𝐤^{})𝐑}\underset{¯}{G}_𝐤^0(\omega )\underset{¯}{T}(\omega )\underset{¯}{G}_𝐤^{}^0(\omega ),`$ (1)
where $`\underset{¯}{T}`$ is the exact $`T`$-matrix for the impurity. We now ask what the change in the current response due to a single impurity is in this system. Note that we intend to study at first only the linear response to an external vector potential $`𝐀`$ à la Kubo, but in the “unperturbed state” described by the exact one-particle eigenstates with the single impurity present. This is simply the usual nonlocal expression for the current response in Coulomb gauge,

$`𝐉(𝐫)={\displaystyle 𝑑𝐫^{}\underset{¯}{K}(𝐫,𝐫^{})𝐀(𝐫^{})},`$ (2)
with $`\underset{¯}{K}(𝐫,𝐫^{})=c/(4\pi \lambda ^2(T))\delta (𝐫𝐫^{})+\underset{¯}{\delta K}(𝐫,𝐫^{})`$. Here we take $`\lambda (T)`$ to be the unperturbed London penetration depth of the $`d`$-wave superconductor, with pure behavior $`\lambda (T)\lambda _0(1+(\mathrm{log}2)T/\mathrm{\Delta }_0)`$ for $`T\mathrm{\Delta }_0`$ and $`\lambda _0\sqrt{mc^2/4\pi ne^2}`$. The Fourier transform $`\underset{¯}{\delta K}(𝐩,𝐪)`$ is now the change in response due to the impurity, which may be easily expressed in terms of the exact 1-impurity $`T`$-matrix using (1),

$`\delta K^{\alpha \beta }(𝐩,𝐪)={\displaystyle \frac{e^2}{c𝒱}}e^{i(𝐩𝐪)𝐑}T{\displaystyle \underset{\omega }{}}{\displaystyle \underset{𝐥}{}}{\displaystyle \frac{1}{2}}\mathrm{Tr}`$ (3)
$`\{𝐯_{𝐥\frac{𝐪}{2}}^\alpha 𝐯_{𝐥\frac{𝐩}{2}}^\beta \underset{¯}{G}_𝐥^0(\omega )\underset{¯}{G}_{𝐥𝐪}^0(\omega )\underset{¯}{T}(\omega )\underset{¯}{G}_{𝐥𝐩}^0(\omega )`$
$`+𝐯_{𝐥+\frac{𝐪}{2}}^\alpha 𝐯_{𝐥+\frac{𝐩}{2}}^\beta \underset{¯}{G}_{𝐥+𝐩}^0(\omega )\underset{¯}{T}(\omega )\underset{¯}{G}_{𝐥+𝐪}^0(\omega )\underset{¯}{G}_𝐥^0(\omega )+{\displaystyle \frac{1}{𝒱}}{\displaystyle \underset{𝐤}{}}`$
$`𝐯_{𝐥+\frac{𝐪}{2}}^\alpha 𝐯_{𝐤\frac{𝐩}{2}}^\beta \underset{¯}{G}_𝐤^0(\omega )\underset{¯}{T}(\omega )\underset{¯}{G}_{𝐥+𝐪}^0(\omega )\underset{¯}{G}_𝐥^0(\omega )\underset{¯}{T}(\omega )\underset{¯}{G}_{𝐤𝐩}^0(\omega )\},`$
where $`\omega `$ is an internal fermion Matsubara frequency, $`𝒱`$ is the volume of the system, and $`𝐯_𝐤ϵ/𝐤`$ is the electron velocity.

We now need to solve (2) together with Maxwell’s equations to determine the spatial dependence of the vector potential caused by the combined effects of Meissner screening and impurity scattering. Since we are primarily interested in long-wavelength effects, the perturbation may be considered to be of order $`1/N`$, where $`N`$ is the number of atoms in the crystal; we may therefore clearly treat the problem perturbatively, by writing $`𝐀(𝐫)=𝐀_0(z)+\delta 𝐀(𝐫)`$, where the unperturbed solution is taken to be the London result (we assume $`\lambda _0>>\xi _0`$ for the homogeneous system) $`𝐀_0(z)=𝐀_0(0)\mathrm{exp}|z|/\lambda _0`$. Note $`z>0`$ is the coordinate describing the theorist’s half-space of superconducting material, and the solution $`A(𝐫)`$ is extended to unphysical values $`z<0`$ as an “image vector potential” to allow a solution by Fourier transform. Specular scattering of quasiparticles from the actual surface is assumed throughout.

To linear order in $`\delta 𝐀`$, the problem can be cast as a differential rather than integro-differential equation,

$`^2\delta 𝐀{\displaystyle \frac{1}{\lambda (T)^2}}\delta 𝐀={\displaystyle \frac{4\pi }{c}}{\displaystyle 𝑑𝐫^{}\underset{¯}{\delta K}(𝐫,𝐫^{})𝐀_0(z^{})}`$ (4)
whose solution can be obtained by Fourier transform,

$`\delta 𝐀(𝐫)`$ $`=`$ $`{\displaystyle \underset{𝐪}{}}e^{i𝐪𝐫}{\displaystyle \frac{𝐅(𝐪)}{q^2+1/\lambda ^2}}`$
$`𝐅(𝐪)`$ $`=`$ $`{\displaystyle \frac{4\pi }{c}}{\displaystyle 𝑑𝐫e^{i𝐪𝐫}𝑑𝐫^{}\underset{¯}{\delta K}(𝐫,𝐫^{})𝐀_0(z^{})}.`$ (5)
With the expressions (2) and (5) it is clearly straightforward, if tedious, to evaluate the full spatial dependence of the vector potential or the current around the impurity site. It is clear from physical considerations that the current at any given point in space $`𝐫`$ will depend both on the distance from the surface $`z`$ and on the distance from the impurity $`|𝐫𝐑|`$. In general, then, the currents and field $`𝐁=\times 𝐀`$ are functions of $`x,y`$ and $`z`$ everywhere. We postpone the full evaluation to a subsequent study. At present we would like simply to obtain an estimate of the size of this effect, which we do by examining the perturbation of the component of the field along the initial applied field $`𝐁(0)`$ which we take to be along $`y`$. The penetration depth is usually defined, even in cases where the decay of the fields is not exponential, i.e. nonlocal electrodynamics, as $`\lambda _{tot}_0^{\mathrm{}}𝑑zB_y/B_y(0)`$. So for our current purposes we will simply adopt this definition and evaluate $`B_y(x=y=0,z)`$, i.e. at the transverse position of the impurity $`𝐑=(0,0,z_0)`$. In this approximation we find

$`{\displaystyle \frac{\delta \lambda }{\lambda }}{\displaystyle \frac{\delta A_x(0)+_0^{\mathrm{}}𝑑z_x\delta A_z+\lambda _z\delta A_x(0)\lambda _x\delta A_z(0)}{A_{0x}(0)}}`$ (6)
A careful examination of these terms and glance at Eq. (5) shows that the size of $`\delta (1/\lambda ^2)2\delta \lambda /\lambda ^3`$ is simply set by $`\underset{¯}{\delta K}(𝐩,𝐪)`$ with $`p,q1/\lambda `$, as intuitively expected. Expanding Eq. (3) for $`p,qk_F`$ and performing the integration over energy $`\xi _𝐤`$ yields from the first two terms

$`\delta K^{xx}(𝐩,𝐩){\displaystyle \frac{ic}{4\pi \lambda _0^2}}T{\displaystyle \underset{\omega }{}}{\displaystyle \frac{T_0}{𝒱}}{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle \frac{2\omega \mathrm{\Delta }_𝐤^2(12\xi _1^2+\eta _𝐩^2)}{\xi _1^3(4\xi _1^2+\eta _𝐩^2)^2}},`$ (7)
where $`\xi _1\sqrt{\omega ^2+\mathrm{\Delta }_𝐤^2}`$, $`T_0=(\pi N_0)^1G_0/(c^2G_0^2)`$, with $`G_0=(1/2\pi N_0)\mathrm{Tr}_𝐤\underset{¯}{G}_𝐤^0=(2i\omega /\pi \xi _1N_0)𝒦(\mathrm{\Delta }_0/\xi _1)`$, is the $`\tau _0`$ component of the $`T`$-matrix, and $`c`$ is the cotangent of the $`s`$-wave scattering phase shift ($`c=0`$ corresponds to infinitely strong scattering). $`𝒦`$ is the complete elliptic integral of the first kind. Note that the third term of (3) contains two factors of $`T_0^2`$ and therefore formally of order $`1/N^2`$; we therefore neglect it in its effect on long-wavelength properties like $`\delta \lambda `$, although it will contribute to local properties.

The quantity $`\eta _𝐩v_𝐤𝐩v_F^z/\lambda `$, where the last approximation follows because the primary spatial gradients are perpendicular to the surface with normal $`\widehat{z}`$, and because we are interested in typical $`p1/\lambda `$. In Figure 1, we show a numerical evaluation of (7) for three different values of the impurity scattering phase shift. As the scattering approaches the unitarity limit, the perturbation to the penetration depth is seen to diverge with decreasing temperature.

To get an analytical estimate of the upturn in the unitarity limit, we first perform a nodal expansion of the order parameter, $`\mathrm{\Delta }_𝐤2\mathrm{\Delta }_0(\varphi \varphi _n)`$, where $`\varphi _n\pi /4`$, $`3\pi /4`$, …, leading for $`T\mathrm{\Delta }_0`$ to the result

$`\delta K^{xx}(𝐩,𝐩)`$ $``$ $`{\displaystyle \frac{c}{4\pi \lambda _0^2}}{\displaystyle \frac{N_0^1}{2\pi ^2T𝒱}}{\displaystyle \frac{1}{n_0\mathrm{log}\frac{2\mathrm{\Delta }_0}{\pi n_0T}}},`$ (8)
where the infrared divergence in (7) was cut off by $`\eta _𝐩`$ or $`c\mathrm{\Delta }_0`$, and $`n_0`$ is the appropriate minimum Matsubara index $`n_0\mathrm{max}(\eta _𝐩/\pi T,c\mathrm{\Delta }_0/\pi T,1)`$. For a (001) surface in the cuprates, $`\eta _𝐩(v_F^z/v_F^{})(\xi _0/\lambda _0)\mathrm{\Delta }_010^310^4`$, so we get a divergence which goes as $`\delta \lambda \lambda _0E_F/[NT\mathrm{log}(\mathrm{\Delta }_0/T)]`$. The factor of $`1/N`$ arises since (8) is proportional to $`N_0^1/𝒱a^3/𝒱`$, where $`a`$ is the lattice spacing. The estimate holds only for temperatures down to $`T_1^{}\mathrm{max}(\eta _𝐩,c\mathrm{\Delta }_0)`$ at which the divergence is cut off, of order tenths of a Kelvin or less in the cuprates if $`c`$ is taken to be 0.

Finite density of scatterers. We would now like to extend the above estimate for a single impurity to the case when a finite density of impurities is present. As mentioned above, when brought into proximity these impurities interfere with one another via hybridization of quasiparticle bound state wave functions, leading eventually to an “impurity band” and residual density of states at the Fermi level, at least in 3Dpjhjltp . In this situation the SCTMA is the appropriate approximation, and a $`T^2`$ behavior in the penetration depth is found rather than an upturn.grossetal  The conditions for the formation of this band are subtle, and we do not address these questions here.

Intuitively, it seems likely that impurities will act as independent sources of current distortion if they are sufficiently far from one another, but it is not so simple to specify what “sufficiently far” means. We note that both the distortion of the current distribution arising from $`\delta K(𝐫,𝐫^{})`$ and the impurity bound state wave functions generically decay over lengths of a few $`\xi _0`$. If the typical interimpurity spacing $`\mathrm{}n_i^{1/2}a`$ (in 2D) is much greater than $`\xi _0`$, each impurity in a penetration length $`\lambda `$ contributes $`\delta \lambda `$ to the observed $`\delta \lambda _{tot}`$. On the other hand, even if $`\mathrm{}<<\xi _0`$, there is a probability $`P(R)\mathrm{exp}R^2/2\mathrm{}^2`$ that an impurity will find itself with no other impurities within a radius $`R`$, if all impurities are randomly distributed. We therefore expect that an upper bound to the observable penetration depth upturn in the case of randomly distributed unitarity limit scatterers will be

$`\delta \lambda _{tot}{\displaystyle \frac{n_iP(\xi _0)E_F}{4\pi ^2\stackrel{~}{T}\mathrm{log}(2\mathrm{\Delta }_0/\pi \stackrel{~}{T})}}\lambda _0,`$ (9)
with $`\stackrel{~}{T}=T+T_1^{}`$. The estimate is an upper bound because we do not currently know how to calculate accurately a geometric factor arising from the contribution of those impurities which are located such that the tails of nodal quasiparticle resonant states (which do not decay as $`e^{\xi _0/r}`$ but rather as $`1/r`$ at $`\omega =0`$)balatsky ; yashenkin  overlap. This factor will lead to significant hybridization and smearing of the singularity at impurity densities below that given by the criterion above. We nevertheless substitute some numbers into (9) to check for plausibility, taking $`E_F/\mathrm{\Delta }_010`$, $`\lambda =1500`$Å, and $`\xi _0=25`$Å for the cuprates. In Figure 2, we plot the results of (9) in laboratory units given these assumptions.

It is seen that, within this crude approximation, the size of the upturn expected is not monotonic with impurity concentration, but at a temperature of 1K peaks around $`n_i10^3`$. The position $`T_2^{}`$ of the minimum in $`\lambda _{tot}(T)`$ will now be set by comparing the impurity-induced change in $`\lambda `$ with that due to thermally excited quasiparticles, $`\delta \lambda (n_i,T_2^{})T_2^{}\lambda _0/\mathrm{\Delta }_0`$; for impurity concentrations of $`n_i=10^310^2`$, this temperature will be of order 1-2K, as seen in Figure 2.

The size of the increases predicted seem quite reasonable, and are of the same order of magnitude as the upturns seen in the work of Bonn et al.bonnetal  and Carrington et al.prozorov  Although the latter authors observed a subtantial increase in the size of the signal when the proportion of (110) surface in the given sample was maximized, indicative of an Andreev surface bound state contribution, there was a substantial signal even in samples with only (100) surfaces, where surface bound states should not exist. Although the authors’ suggestion that this result is due to (110) faceting in these samples is possible, the mechanism we suggest may also be present.

While the estimate (9) is based on the assumption of randomly distributed impurities, it is important to recognize that in real systems clustering will take place. In this case even a relatively disordered system may have significant numbers of isolated impurities or isolated atomic-scale clusters which give rise to a singular current response. In this case we might expect significant sample-to-sample variation in $`\delta \lambda `$ in samples with identical average concentration $`n_i`$. This appears to be precisely the effect found by Bonn et al. in their measurements on YBCO single crystals with 0.31% Znbonnetal .

Magnetic field. Since the penetration depth upturn in the case considered here is due to the large number of quasiparticle excitations near the nodal directions, as is the similar upturn in the case of Andreev surface states, we might expect any physical effect which smears the gap nodes to cut off the upturn. In particular, the orbital coupling to an applied magnetic field (nonlinear electrodynamics) will suppress the upturn as it does in the Andreev casebarashetal . We need only add the Doppler shift of the quasiparticle energy in any of the expressions above, $`i\omega i\omega +𝐯_s𝐤`$, where $`𝐯_s`$ is the local superfluid velocity. A typical shift $`v_sk_F(H/H_0)\mathrm{\Delta }_0`$, where $`H`$ is the applied field and $`H_0=3\mathrm{\Phi }_0/(\pi ^2\xi _0\lambda _0)`$ is of the order of the thermodynamic critical field, and the effect of the field may be included approximately by generalizing $`\stackrel{~}{T}T+T_1^{}+(H/H_0)\mathrm{\Delta }_0`$ in (9). We postpone a more quantitative study of the field dependence to a later work.

Conclusions. We predict that a $`d`$-wave superconductor with isolated strong nonmagnetic impurities should exhibit an upturn in the penetration depth varying as $`1/T\mathrm{log}T`$ below a temperature scale set by the disorder in the system and above one set by the bulk penetration depth or the scattering phase shift. This effect may have contributed to upturns already observed in experiments where Andreev surface states, which produce a similar effect, are not indicatedbonnetal ; prozorov . The physical origin of the upturn in both cases is the depletion of the supercurrent at low $`T`$ in low-energy bound states, but the impurity-induced upturn is independent of surface geometry. We argue that the magnitude of the effect should have a characteristic dependence not only on the concentration of impurities, but on the exact nature of the statistical distribution of impurities in the sample. We emphasize that the estimates presented here are merely a crude plausibility argument that if such upturns are observed at low temperatures, they could arise from singular current response of local defects in the crystal, and that it will be useful to study the dependence on disorder. Our work also allows in principle an exact calculation of the local current flow around an impurity; it will be interesting to see how the singular quasiparticle currents are distributed in the neighborhood of the impurity site, and whether these features can be detected by STM.

Acknowledgements. This work is supported by NSF DMR-9974396, by BMBF 13N6918/1, and by a grant from the A. v. Humboldt foundation. The authors are grateful to Jochen Mannhart and the EP6 Lehrstuhl at the University of Augsburg for hospitality during the completion of this work, and to Russ Giannetta for a valuable communication.


format: text

text:
# On directed interacting animals and directed percolation

## 1 Introduction

Ever since the seminal work of Broadbent and Hammersley in 1957 , the lasting interest in directed percolation and related statistical models has been nurtured by several motivations. Firstly, very anisotropic fractal structures are observed in Nature in a variety of situations where a bias is introduced by an external field or a gradient, ranging from river basins  to stress cracks in metals  and to lightning discharges. The preferred direction may alternatively correspond to time, with fractal patterns appearing in the dynamics of the systems considered. Pomeau  thus proposed that in simple models of turbulence the spatio-temporal structure of the active regions could be described by directed percolation (DP in the following), which was later confirmed by extensive numerical simulations . More generally the question arises if the large-scale properties of the observed structures depend on the details of the dynamical processes occurring on short scales, or if they are more universal and can be described in terms of purely statistical models.

Another motivation is that directed systems may be studied using renormalization group (RG) tools, providing a possible approach to nonequilibrium phenomena. The introduction of a preferred direction is a relevant perturbation in the RG sense, it is expected to change the universality class of a phase transition and modify the critical exponents in a non-trivial way. In particular the DP threshold is in the same universality class as Reggeon field theory . In spite of much work, however, exact results are known only for simpler variants and DP itself remains a challenge for theorists, being one of the few fundamental problems to remain unsolved in two dimensions \- see  for a recent review of this model and its position among models of non-equilibrium phase transitions. As pointed out by Grassberger  it is surprising in view of the central role played by DP that no detailed experimental tests of the predictions for its critical behaviour have been carried out. A serious difficulty is that in practice many effects can obscure the expected behaviour at attainable scales, so it is of interest to explore the crossover regime between DP and other universality classes.

Directed lattice animals provide another class of fractal structures that may be relevant for oriented systems, and it has been suggested for example that they might describe the large-scale geometry of river networks . The basic animal model consists in enumerating all connected clusters of N sites on a lattice but various extensions may be considered by introducing local interactions between the cluster sites. One can in this way obtain different statistical ensembles ranging from loopless structures (lattice trees) to fully compact ones. In contrast with the situation for DP several exact results are known for directed animals on various lattices  \- for recent results and a review of what is known rigorously see  and references therein.

A remarkable connection between the two types of systems has been found by Dhar , who showed that a model of interacting directed animals (IDA) can be mapped onto a directed percolation problem if certain relations hold between their respective parameters. For sufficiently strong attractive interactions the typical animals correspond to percolation clusters in the dense phase above the percolation threshold. As a consequence the collapse transition of the IDA model is expected to be of the same nature as the DP threshold, in contrast to the isotropic case where the percolation threshold has been found to correspond to a higher-order multi-critical point in the phase diagram of interacting animals .

We consider in the present work the statistical properties of fully directed lattice animals with nearest-neighbour interactions on the square lattice, using a transfer matrix approach to calculate their correlation lengths on infinite strips of finite width. The results are analysed using a phenomenological renormalization group method adapted to directed systems , which allows an accurate determination of the critical fugacity.

The paper is organized as follows: in the next section the IDA model and the main quantities of interest are introduced and we present the correspondences with a percolation problem and with the generating functions of various simple restricted ensembles. In section 3 we apply the phenomenological renormalization (PR) group approach to non-interacting animals and check its accuracy by comparing the results obtained using strips of widths up to $`n=18`$ with the exact values or the best available ones. Section 4 is devoted to the study of the behaviour at the collapse transition temperature, we determine the corresponding critical geometrical and thermal exponents and compare them with DP exponents. The low-temperature regime, which corresponds to the percolating phase of DP, is considered in section 5. In this regime the transition occurring at the critical fugacity is first-order and the grand partition function displays an essential singularity, but the finite-size PR analysis remains applicable and yields results in excellent agreement with Dhar’s prediction for the critical line . In section 6 we obtain very precise values for the critical fugacities of various ensembles, such as bond animals, site or bond trees and strongly embeddable animals , and we compare them with published values when these exist. In section 7 particular attention is paid to the crossover region in the immediate vicinity of the collapse transition. Though the singularity is very weak we are able to obtain a direct estimate of the crossover exponent. We present analytic arguments which show that this crossover exponent is equal to the Fisher exponent $`\sigma `$ for DP , in agreement with our numerical results.

## 2 The interacting animal model

### 2.1 Directed animals

A directed animal is a connected cluster of $`N`$ sites, some of which belong to the root, such that any other site of the animal can be reached from the root by a walk which never goes opposite to the directed axis (see figure 1). In the simplest case, non-interacting animals, one associates the same weight $`x^N`$ to each cluster having $`N`$ sites, so the basic generating function is

$$S_0(x)=\underset{N}{}\mathrm{\Omega }_Nx^N,$$
(1)
where $`\mathrm{\Omega }_N`$ is the number of directed connected clusters having a single site (the origin) as root.

The problem of enumerating such objects on lattices was recognized as an interesting combinatorial problem and received much attention over the last two decades, in part because it turned out to be possible to obtain many exact results concerning the $`\mathrm{\Omega }_N`$, whereas for undirected animals only bounds are known. In particular the generating function (1) is known for the square lattice studied here :

$$S_0(x)=\frac{1}{2}((\frac{1+x}{13x})^{1/2}1).$$
(2)
For physical applications the metric properties of directed animals are perhaps yet more important. They are best described in terms of two independent characteristic lengths with asymptotic behaviour of the form

$$R^{}c_1N^\nu ^{},R^{}c_2N^\nu ^{},$$
(3)
where $`R^{}`$ provides a measure of the average distance between the root and the most distant sites along the directed axis, while $`R^{}`$ corresponds to the average distance between two extreme animal sites in the perpendicular (transverse) direction. The value of the perpendicular exponent is known to be $`\nu ^{}=1/2`$ in $`d=2`$, from the exact solution for various lattices  and in agreement with a field theory argument connecting the large-scale behavior of directed animals and the dynamics at the Yang-Lee edge singularity . The longitudinal exponent $`\nu ^{}`$ has proven much more difficult to obtain and is known only numerically .

### 2.2 Branched polymers and interacting animals

Interacting animals are often used as simple lattice models to describe the behaviour of branched polymers in dilute solutions. Indeed, field theory calculations  indicate that the existence of loops in animals is irrelevant in the RG sense, so one expects them to belong to the same universality class as branched polymers and to share the same large-scale behaviour. A simple way to take into account the interactions with the solvent is to associate an energy $`ϵ`$ with every pair of neighbouring occupied sites. This energy corresponds to the difference between the monomer-solvent and the monomer-monomer interaction energies (one can distinguish between chemically bonded pairs and others in mere physical contact , but the present simple model already contains much of the physics). At an absolute temperature $`T`$ the Boltzmann thermal weight associated with each pair is then

$$w=\mathrm{exp}(ϵ/k_BT),$$
(4)
where $`k_B`$ is the Boltzmann constant.

A partition function for polymers made of exactly $`N`$ monomers may be defined:

$$Z_N(T)=\underset{P}{}\mathrm{\Omega }(N,P)w^P,$$
(5)
where $`\mathrm{\Omega }(N,P)`$ is the number of clusters of $`N`$ sites having $`P`$ pairs of occupied nearest-neighbors (note that all possible polymer topologies contribute to $`Z_N`$, this corresponds to a model of ”living polymers” where the monomers may rearrange). At infinite temperature $`w=1`$ and one recovers the statistics of an ensemble of equally weighted animals. A good solvent corresponds to $`ϵ0`$ ($`w1`$), a poor solvent corresponds to $`ϵ>>0`$ (large $`w`$). The free energy per site in the thermodynamic limit is given by

$$f(T)=\underset{N\mathrm{}}{lim}\frac{k_BT}{N}\mathrm{log}Z_N(T).$$
(6)
It is often convenient to consider the generating function

$$G(x,w)=\underset{N}{}Z_N(w)x^N=\underset{NP}{}\mathrm{\Omega }(N,P)x^Nw^P,$$
(7)
which is analogous to a grand partition function and has a singularity for $`x=x_c(w)`$, the critical fugacity. From (5), (6) and (7) one has

$$f(T)=k_BT\mathrm{log}x_c(w).$$
(8)
### 2.3 Collapse transition and percolation threshold

When the temperature is lowered compact configurations become more probable and the polymers are expected to shrink from an extended state to a dense globular state. In the limit $`N\mathrm{}`$ a collapse transition occurs between these two states at a well-defined temperature (for a statistically typical polymer). This transition is described by a critical point, usually referred to as the $`\theta `$-point, which is analogous to a tricritical point for magnetic systems. It is believed that the average linear size of the polymers at the $`\theta `$-point scales as $`N^{\nu _\theta }`$, where $`\nu _\theta `$ differs from the exponents $`\nu `$ and $`\nu _c`$ corresponding respectively to the extended and dense states. Various lattice models have been proposed for the collapse of isotropic branched polymers . No complete solution of these models on regular lattices is known, and most of them have been studied by numerical methods (see, e.g., ) such as Monte Carlo simulation, exact enumeration on lattices, or transfer-matrix calculations coupled with phenomenological renormalization (PR) .

For directed animals the corresponding transition was first studied by Monte Carlo simulations , within an isotropic scaling approach. Although this assumption did not properly take into account the intrinsically anisotropic nature of directed animals, the value obtained for the $`\theta `$ temperature is found to be in agreement with more accurate approaches.

A central result, due to Dhar , and rederived independently in a recent study of the adsorption of interacting animals , is the existence of a correspondence between the partition function of directed lattice animals with interactions between first and second neighbours and a site-bond percolation problem on the same directed square lattice. In particular the IDA first-neighbour model studied here corresponds to a simple bond percolation problem: a given animal $`A`$ occurs with probability

$$𝒫(A)=\frac{1+q}{p}(\frac{pq^2}{1+q})^N(\frac{1+q}{q})^P$$
(9)
and using (5) one obtains the relation

$$Pr(N)=\underset{|A|=N}{}𝒫(A)=\frac{1+q}{p}(\frac{pq^2}{1+q})^NZ_N(w)$$
(10)
where the sum bears over all clusters of $`N`$ sites and

$$p=1q=(w2)/(w1)$$
(11)
is the associated bond percolation probability. $`Pr(N)`$ is the probability that for this value of $`p`$ a percolation cluster starting from the origin contains precisely $`N`$ sites. Hence in the thermodynamic limit the free energy per site is given for $`w>2`$ by

$$\frac{f(w)}{k_BT}=\mathrm{log}(\frac{pq^2}{1+q})\underset{N\mathrm{}}{lim}(1/N)\mathrm{log}Pr(N),$$
(12)
For isotropic percolation in $`d`$ dimensions it has been proved  that above the percolation threshold the decay of $`Pr(N)`$ is slower than exponential and that there exists a value $`p_1>p_c`$ such that for $`p>p_1`$

$$\mathrm{log}Pr(N)CN^{11/d}\text{if}N\mathrm{},$$
(13)
but in fact this behaviour is expected to hold for all $`p>p_c`$. Physically this just means that most sites of a large finite cluster lie far from its external boundary and belong to its ”bulk”, whose effective volume is of the order of the geometric one, while the external boundary that needs to be disconnected from the infinite cluster is not fractal and contains of the order of $`N^{(d1)/d}`$ sites .

In the directed case large finite percolation clusters are quite rare above $`p_c`$ and difficult to study directly, but a proof has been outlined recently that in this case too the decay of $`Pr(N)`$ is slower than exponential  and it is natural to expect that $`\mathrm{log}Pr(N)`$ also follows a power law decay for all $`p>p_c`$. The second term on the right-hand side of (12) is then negligible in the thermodynamic limit and one obtains a simple equation for the critical line in the low-temperature phase :

$$x_c(w)=\frac{pq^2}{1+q}=\frac{w2}{w(w1)^2}\text{if}w>w(\theta ).$$
(14)
A remarkable consequence is that the free energy per site is not singular on the low-T side of the transition. The excellent agreement of this relation with our numerical results presented below confirms its validity. As a further consequence the critical temperature of the IDA model studied here is related to the critical probability $`p_c`$ for the directed bond percolation process on the same lattice:

$$w(\theta )=(2p_c)/(1p_c).$$
(15)
The threshold value $`p_c`$ as well as the values of DP critical exponents have recently been estimated very accurately using low-density series expansions , so the thermal critical exponents of collapsing directed animals are also known with high precision.

Note, however, that the IDA model is more general than DP as it remains well defined for repulsive interactions ($`w<1`$), whereas the equivalence with DP is only valid for $`w>2`$. Also, the whole $`(x,w)`$ plane may be given a physical meaning if the region $`x>x_c(w)`$ is interpreted as a ”weak gel” stabilized through application of an external pressure .

### 2.4 Special ensembles

Let us remark at this stage that the model defined by (7) contains in its phase diagram several special points related to ensembles defined by simple geometrical properties.

Note firstly that for directed animals the number $``$ of independent loops is just equal to the number $`N_2`$ of sites connected to two occupied predecessors. Every such site contributes two pairs of occupied neighbours in (7) while the other occupied sites contribute only one, except for the origin (for simplicity we consider animals rooted at a single site), so the total number of pairs is $`P=2N_2+(N1N_2)`$. The number of loops is then

$$=PN+1$$
(16)
and the loop generating function  is given by

$`G_l(y,v)`$ $`=`$ $`{\displaystyle \underset{N}{}}\mathrm{\Omega }(N,)y^Nv^{}`$ (17)
$`=`$ $`vG(y/v,v).`$ (18)
In particular one obtains the generating function $`G_{st}(y)`$ of *site trees* (loopless animals):

$$G_{st}(y)=G_l(y,0)=\underset{w0,x\mathrm{}}{lim}wG(x,w)$$
(19)
with $`y=wx`$ fixed.

A simple correspondence also exists between non-interacting directed *bond animals* (i.e. clusters of connected directed bonds) and site animals: For every occupied site having only one predecessor consider the bond linking them, while for the $`N_2`$ sites with two predecessors consider one of three possibilities, either chosing one of the connecting bonds or keeping both of them. Since $`N_2`$ is equal to the number $``$ of loops this correspondence implies that the generating function of these bond animals is given by

$`G_b(y)`$ $`=`$ $`{\displaystyle \underset{N}{}}\mathrm{\Omega }(N,)y^{N1}(2y+y^2)^{}`$ (20)
$`=`$ $`{\displaystyle \frac{1}{y}}G_l(y,2+y)`$ (21)
$`=`$ $`{\displaystyle \frac{1}{x}}G(x,{\displaystyle \frac{2}{1x}}),`$ (22)
with $`x=y/(2+y)`$ and where we have used (17) and (18) to obtain (21) and (22) respectively.

If one now discards from the previous bond configurations all the loop-forming ones (i.e. those with two connecting bonds going to the same site) one is left with directed *bond trees*. The correspondence with site animals just discussed shows that their generating function is given by

$`G_{bt}(y)`$ $`=`$ $`{\displaystyle \frac{1}{y}}G_l(y,2)`$ (23)
$`=`$ $`{\displaystyle \frac{1}{x}}G(x,2),`$ (24)
with $`x=y/2`$.

Finally, *strongly embeddable* animals  correspond to maximally connected bond configurations, in the sense that for a given site animal one systematically keeps the two bonds going to every site with two predecessors (bond trees correspond to minimally connected configurations). Their generating function is then simply

$`G_{emb}(y)`$ $`=`$ $`{\displaystyle \underset{N}{}}\mathrm{\Omega }(N,)y^{N1}(y^2)^{}`$ (25)
$`=`$ $`{\displaystyle \frac{1}{y}}G_l(y,y)`$ (26)
$`=`$ $`G(1,y).`$ (27)
Detailed numerical results for these various ensembles are presented in section 6.

## 3 Numerical methods and phenomenological renormalization analysis

### 3.1 Transfer matrix method

Here we shall use the method of transfer matrices coupled to a phenomenological renormalization analysis , that was applied to lattice animals in . It is convenient to define the generating function $`G_{0𝐑}(x,w)`$ for the total number $`\omega _{0𝐑}(N,P)`$ of directed animals of $`N`$ sites with $`P`$ pairs of nearest-neighbors connecting the points $`0`$ and $`𝐑`$

$$G_{0𝐑}(x,w)=\underset{NP}{}\omega _{0𝐑}(N,P)x^Nw^P.$$
(28)
On a strip of width $`n`$ with, say, periodic boundary conditions, linear recursion relations may be written between the $`\omega _{0𝐑}(N,P)`$ for two successive values of the distance $`R`$ along the strip. The asymptotic behaviour of $`G_{0𝐑}(x,w)`$ is dominated by the largest eigenvalue $`\lambda _n(x,w)`$ of the corresponding transfer matrix $`𝒯_n`$ :

$$G_{0𝐑}(x,w)[\lambda _n(x,w)]^R=\mathrm{exp}\left(\frac{R}{\xi _n^{}(x,w)}\right),R\mathrm{},$$
(29)
which defines the longitudinal (or parallel) correlation length:

$$\xi _n^{}(x,w)=1/\mathrm{log}\lambda _n.$$
(30)
This correlation length diverges along the critical line $`x=\stackrel{~}{x}_n(w)`$, where $`\stackrel{~}{x}_n`$ is the smallest positive root of

$$\lambda _n(\stackrel{~}{x}_n(w),w)=1$$
(31)
and the generating function (28) is infinite for $`x>\stackrel{~}{x}_n(w)`$. Let us recall that this line determines the free energy per site of a very large animal of width $`n`$, $`f_n=k_BT\mathrm{log}\stackrel{~}{x}_n(w)`$. In a similar way $`\lambda _n(x,w)`$ determines the average site density $`\rho _n(w)`$ of an animal on the strip  through

$$\rho _n(w)=\frac{1}{n}\frac{\mathrm{log}[\lambda _n(x,w)]}{\mathrm{log}x}|_{x=\stackrel{~}{x}_n},$$
(32)
where the derivative is calculated at $`\stackrel{~}{x}_n(w)`$.

The gel phase, for $`x>\stackrel{~}{x}_n(w)`$, will not be studied in the present work, let us just note that expression (32) remains valid and gives the gel density , but the expression of the correlation length is different (see  for a lucid discussion of that point).

### 3.2 Phenomenological renormalization in the high - temperature limit

In the high-temperature limit ($`w=1`$) one recovers non-interacting lattice animals which as mentioned above have been studied by a variety of methods (see  and references cited therein). Let us recall here just two recent results. Using a connection of this problem to the continuous-time dynamics of a one-dimensional lattice gas , Dhar was able to obtain a very accurate value of the longitudinal size exponent $`\nu ^{}=\mathrm{0.817\hspace{0.17em}33}(5)`$, where $`(5)`$ indicates the estimated uncertainty for the last digit. This critical exponent was estimated more recently  by the method of series expansions, which gave $`\nu ^{}=\mathrm{0.817\hspace{0.17em}22}(5)`$. Both these results exclude the conjecture $`\nu ^{}=9/11`$ which was based on phenomenological renormalization on strips of width up to $`n=12`$ , and it is now generally believed that $`\nu ^{}`$ is not a rational number.

We shall first extend the latter calculation on strips of width up to $`n=18`$. As we shall show, our estimate of $`\nu ^{}`$ is in excellent agreement with the results just quoted, which provides a useful check of the overall reliability of the present approach to the case of interacting animals.

According to renormalization-group theory it is expected that for $`n1`$ and $`xx_c1`$ the longitudinal correlation length satisfies a finite-size scaling relation of the form

$$\xi _n^{}n^\kappa F[n^{1/\nu ^{}}(xx_c)],$$
(33)
where the scaling exponent $`\kappa =\nu ^{}/\nu ^{}`$ and the scaling function $`F(z)`$ is finite and regular for $`z=0`$. This allows to obtain an estimate of the critical fugacity $`x_c`$ as the fixed point of phenomenological renormalization equations involving three different strips. Using three consecutive widths $`(n1,n,n+1)`$ the equations for the fixed point read

$$\frac{\xi _{n1}^{}(x_{c,n})}{(n1)^{\kappa _n}}=\frac{\xi _n^{}(x_{c,n})}{n^{\kappa _n}}=\frac{\xi _{n+1}^{}(x_{c,n})}{(n+1)^{\kappa _n}},$$
(34)
where $`x_{c,n}`$ and $`\kappa _n`$ are respectively estimates of the critical fugacity and the scaling exponent. Eliminating $`\kappa _n`$ in (34) gives an equation for $`x_{c,n}`$:

$$\mathrm{log}\frac{\xi _{n+1}^{}}{\xi _n^{}}=\frac{\mathrm{log}\frac{n+1}{n}}{\mathrm{log}\frac{n}{n1}}\mathrm{log}\frac{\xi _n^{}}{\xi _{n1}^{}}\text{for}x=x_{c,n},$$
(35)
which would yield the exact value of $`x_c`$ if the scaling relation (33) held rigorously. Finite-size estimates $`\kappa _n,\nu _n^{}`$ and $`\nu _n^{}`$ of the critical exponents are then obtained from

$`\kappa _n`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}(\xi _{n+1}^{}/\xi _n^{})}{\mathrm{log}[(n+1)/n]}},`$ (36)
$`{\displaystyle \frac{1}{\nu _n^{}}}`$ $`=`$ $`\kappa _n+{\displaystyle \frac{\mathrm{log}[(d\xi _{n+1}^{}/dx)/(d\xi _n^{}/dx)]}{\mathrm{log}[(n+1)/n]}},`$ (37)
$`\nu _n^{}`$ $`=`$ $`\kappa _n\nu _n^{},`$ (38)
where all ratios and derivatives are evaluated at $`x=x_{c,n}`$.

To get accurate numerical estimates of the required derivatives we have used here a perturbative expansion method for the transfer matrices as described in . This enabled us to compute the largest eigenvalue derivatives with an accuracy typically better than $`10^{12}`$. We also took advantage of the fact that the transfer matrix is singular for $`w=1`$ (i.e., $`det(𝒯_n)=0`$) to reduce its size: for instance, for $`n=15`$ the matrix size could be reduced from ($`1223\times 1223`$) to ($`198\times 198`$).

### 3.3 Analysis of the numerical results

In order to extrapolate these finite-size estimates two different extrapolation procedures were used. The first one (to be referred to as extrapolation 1) is based on a simple sequential fit of the form

$$\nu _n^{}=\nu ^{}+An^{\omega _1}+Bn^{\omega _2},$$
(39)
where the estimates of $`\nu ^{}`$ and of $`A,B,\omega _1,\omega _2`$ can be obtained from the values of $`\nu _n^{}`$ for five consecutive widths. One could expect that for large $`n`$ these estimates take certain limiting values, providing the form (39) is well chosen. It turns out, however, that in some extrapolation procedures considered in this work the effective values of $`B`$ and $`\omega _2`$ are rather unstable, which probably indicates that a more complicated form should be used for the next to leading correction. In such cases we prefer to keep only the first two terms on the right hand side of (39), which always seem to lead to a reliable extrapolation. In either case the new sequences of estimates that we calculate are monotonic and they depend on $`n`$ more weakly than the original ones. The final estimates are obtained by plotting the extrapolated data versus $`n^{\stackrel{~}{\omega }}`$ were $`\stackrel{~}{\omega }`$ is adjusted to give a smooth convergence.

The second procedure used to extrapolate for $`n\mathrm{}`$, which will be referred to here as extrapolation 2, is known as the Bulirsch-Stoer (BST) algorithm . This method includes a free parameter $`\omega `$ which has to be adjusted self-consistently to minimize the fluctuations in the sequences of estimates produced by the BST algorithm. Let us note here that the ranges of $`\omega `$ estimated by this method are in agreement with the corresponding estimates for the leading correction to scaling exponent $`\omega _1`$ from (39).

Our finite-size estimates together with the extrapolated values for $`\nu ^{}`$, $`\nu ^{}`$ and $`x_c`$ are presented in table 1. One can notice the excellent agreement of these extrapolated values with the exact results $`\nu ^{}=1/2`$, $`x_c=1/3`$, while the present result $`\nu ^{}=\mathrm{0.817\hspace{0.17em}31}(5)`$ corroborates the best earlier estimates of the longitudinal size exponent.

## 4 $`\theta `$-Point behaviour

### 4.1 Determination of the $`\theta `$ critical temperature

It was shown in  that for non-interacting directed animals the dependence of the average site density (32) on the strip width is of the form

$$\rho _n(w=1)Cn^{\frac{1\nu ^{}}{\nu ^{}}1}$$
(40)
for large $`n`$ (here and in the following we use $`C`$ without super- or underscript as generic notation for an undetermined constant). By analogy with this behaviour and with the case of isotropic lattice animals  this density is expected to have the following finite-size scaling form in the vicinity of the IDA $`\theta `$-point:

$$\rho _n(w)n^{\psi 1}H[(ww(\theta ))n^{1/\nu _t^{}}],$$
(41)
for $`n1`$ and $`ww(\theta )1`$, with

$$\psi =(1\nu _\theta ^{})/\nu _\theta ^{},$$
(42)
where $`\nu _\theta ^{}`$ and $`\nu _\theta ^{}`$ denote respectively the longitudinal and perpendicular size critical exponents at the $`\theta `$-point. $`H(z)`$ is a non-singular scaling function and $`\nu _t^{}`$, the perpendicular thermal correlation length exponent, is expected to be the same as the perpendicular exponent $`\nu _{DP}^{}`$ for directed percolation. Using (32) to calculate $`\rho _n`$ we have determined the quantities

$$\psi _n(w)=\frac{\mathrm{log}(\rho _n/\rho _{n+1})}{\mathrm{log}[n/(n+1)]}+1.$$
(43)
Due to the scaling form (41) one expects that $`\psi _n(w(\theta ))\psi `$ for $`n1`$, so the solutions of the fixed-point equations

$$\psi _{n+1}(w)=\psi _n(w),$$
(44)
yield a sequence of estimates $`w_n(\theta )`$ and $`\psi _{\theta ,n}`$ for the critical weight $`w(\theta )`$ and the ratio $`\psi `$.

The results are given in table 2. Note that this two-parameter phenomenological renormalization (PR) group method gives a rather accurate estimate for the Boltzmann critical weight, $`w(\theta )=\mathrm{3.814\hspace{0.17em}53}(1)`$. This estimate is even better than the one obtained using relation (15) and the result of a one-parameter PR approach for the directed percolation threshold . It is, however, less accurate than the best one: $`w(\theta )=\mathrm{3.814\hspace{0.17em}524\hspace{0.17em}4}`$, deduced from (15) and the value $`p_c=\mathrm{0.644\hspace{0.17em}700\hspace{0.17em}185}(5)`$ recently determined using a one-parameter series expansion method for DP on the square lattice . So, we shall use this last value of $`w(\theta )`$ in order to estimate $`\nu _\theta ^{}`$ and $`\nu _\theta ^{}`$.

### 4.2 Geometrical critical exponents $`\nu _\theta ^{}`$ and $`\nu _\theta ^{}`$

There are different methods to calculate finite-size estimates of two independent critical exponents $`\nu _\theta ^{}`$ and $`\nu _\theta ^{}`$. The simplest one is to suppose that the finite-size scaling form (33) holds for each fixed $`w`$, which allows one to apply the approach described in the preceding subsection. Thus, using three-strip PR for $`w=w(\theta )`$ we have obtained two sequences of estimates $`\nu _{\theta ,n}^{}`$ and $`\nu _{\theta ,n}^{}`$. The results are shown in table 2.

A more sophisticated analysis  consists in taking into account explicitly in the scaling function the two independent relevant fields that control the tricritical behaviour near the $`\theta `$-point. One then writes

$$\xi _n^{}n^{\kappa _\theta }F[n^{1/\nu _1}u,n^{1/\nu _2}v],$$
(45)
where to leading order the scaling fields are related to the model parameters by

$`u`$ $`=`$ $`a(xx_\theta )+b(ww(\theta ))+\mathrm{},`$ (46)
$`v`$ $`=`$ $`c(xx_\theta )+d(ww(\theta ))+\mathrm{},`$ (47)
and

$$\kappa _\theta =\nu _\theta ^{}/\nu _\theta ^{}.$$
(48)
The smaller scaling exponent, $`\nu _1`$, say, is related to the transverse size of the animals at the $`\theta `$-point , so $`\nu _1=\nu _\theta ^{}`$, while the larger one controls the thermal correlation length, so $`\nu _2=\nu _t^{}`$.

Now fixing $`w=w(\theta )`$ one can still use (35) and (36) to obtain a sequence of estimates $`x_{\theta ,n}`$ and $`\kappa _{\theta ,n}`$ for $`x_\theta `$ and $`\kappa _\theta `$ respectively. The extrapolated values for $`x_\theta `$ (see table 3) are in excellent agreement with the value deduced from (14), $`x_\theta =\mathrm{0.060\hspace{0.17em}049\hspace{0.17em}912\hspace{0.17em}4}(10)`$. The main uncertainty comes in fact from the value of $`p_c`$ used to estimate $`w(\theta )`$.

Taking the partial derivatives with respect to $`x`$ and $`w`$ for a given width $`n`$ (we follow closely the method and notations of ), one has at the $`\theta `$-point

$`_x\xi _n^{}`$ $`=`$ $`n^{\kappa _\theta }(an^{1/\nu _1}F/u+cn^{1/\nu _2}F/v),`$ (49)
$`_w\xi _n^{}`$ $`=`$ $`n^{\kappa _\theta }(bn^{1/\nu _1}F/u+dn^{1/\nu _2}F/v).`$ (50)
The unknown constants $`aF/u`$ and $`bF/u`$ may be eliminated using another width $`m`$:

$`{\displaystyle \frac{_x\xi _n^{}}{n^{\kappa _\theta +1/\nu _1}}}{\displaystyle \frac{_x\xi _m^{}}{m^{\kappa _\theta +1/\nu _1}}}`$ $`=`$ $`c(n^{1/\nu _21/\nu _1}m^{1/\nu _21/\nu _1})F/v,`$ (51)
$`{\displaystyle \frac{_w\xi _n^{}}{n^{\kappa _\theta +1/\nu _1}}}{\displaystyle \frac{_w\xi _m^{}}{m^{\kappa _\theta +1/\nu _1}}}`$ $`=`$ $`d(n^{1/\nu _21/\nu _1}m^{1/\nu _21/\nu _1})F/v.`$ (52)
The ratio $`c/d`$ may in turn be eliminated using a third width $`l`$, yielding finally an equation for the exponent $`z=\kappa _\theta +1/\nu _1`$ :

$$\frac{_x\xi _n^{}(n/m)^z_x\xi _n^{}}{_x\xi _n^{}(n/l)^z_x\xi _l^{}}=\frac{_w\xi _n^{}(n/m)^z_w\xi _m^{}}{_w\xi _n^{}(n/l)^z_w\xi _l^{}}.$$
(53)
The same calculation may be performed again exchanging the terms involving $`\nu _1`$ and $`\nu _2`$. One finds that $`z^{}=\kappa _\theta +1/\nu _2`$ also satisfies (53), so one expects this equation to have two positive solutions, from which estimates for $`\nu _1`$ and $`\nu _2`$ may be obtained since $`\kappa _\theta `$ has been estimated independently. The numerical results obtained using three consecutive widths are presented in table 3, where it is readily seen that the extrapolated value $`\nu _t^{}=1.097(1)`$ agrees with the value of the perpendicular percolation exponent $`\nu _{DP}^{}=\mathrm{1.096\hspace{0.17em}854}(4)`$ . The ratio $`\kappa _\theta `$ is also very close to its directed percolation counterpart $`\nu _{DP}^{}/\nu _{DP}^{}=\mathrm{1.580\hspace{0.33em}745}(10)`$ .

This suggests that the correspondence with critical directed percolation clusters is also valid for the geometric properties of IDA clusters in addition to their thermal critical behaviour. Indeed, one finds that $`\nu _\theta ^{}/\nu _{DP}^{}=\mathrm{0.391\hspace{0.33em}55}(7)`$ agrees with the value of the exponent $`\sigma =1/(\beta _{DP}+\gamma _{DP})=\mathrm{0.391\hspace{0.17em}510}(2)`$ . This is in direct analogy with the relation $`\nu _c=\sigma _p\nu _p`$ for isotropic percolation, where $`\nu _c`$ is the exponent associated with the typical radius of large clusters . Although interacting animals are related to directed bond percolation through (10) it was not abvious a priori that their geometrical exponents should also be related simply to those of DP, as they depend on the behaviour of $`G_{0R}(x,w)`$ at fixed $`w`$ and not along Dhar’s line.

### 4.3 Density exponents

Defining now the thermal density exponent $`\beta _t`$ through the variation of the density along the critical line :

$$\rho (T)C(ww(\theta ))^{\beta _t},$$
(54)
one obtains from the scaling relation (41) and the numerical results above for $`\psi `$ and $`\nu _t^{}`$

$$\beta _t=(1\psi )\nu _t^{}=0.2765(4),$$
(55)
in agreement with the percolation exponent $`\beta _{DP}=\mathrm{0.276\hspace{0.33em}486}(8)`$ .

At a fixed temperature, on the other hand, the density in the gel phase (for $`x>x_c(w)`$) is expected to vary close to the critical line as

$$\rho (x)C(xx_c(w))^{\beta _x},$$
(56)
where the exponent $`\beta _x`$ depends on the temperature and is related through scaling relations to the other critical exponents, yielding

$`\beta _x`$ $`=`$ $`\nu _\theta ^{}+\nu _\theta ^{}1=0.1082(1)\text{for}w=w_\theta ,`$ (57)
$`=`$ $`\nu ^{}+\nu ^{}1=0.3173(1)\text{for}w<w_\theta .`$ (58)
We have not checked these predictions numerically.

## 5 Low-temperature region

### 5.1 First-order transition and essential singularity

For $`w>w(\theta )`$ typical animals are in a collapsed state and cover a finite fraction of the lattice sites in the thermodynamic limit. One expects the transition at fixed temperature to be first-order  and the average site density $`\rho (x,w)`$ to jump from zero to a finite value on the critical line given by $`x_c=(w2)/w(w1)^2`$ in that region.

Using relation (10) the animal generating function may be written

$$G(x,w)=\underset{N}{}Z_N(w)x^N=(\frac{w2}{w})\underset{N}{}(x/x_c)^NPr(N).$$
(59)
If we now assume, as discussed in section (2.3), that for all $`p>p_c`$ the probability of large finite directed percolation clusters varies as

$$\mathrm{log}Pr(N)CN^\zeta ,$$
(60)
where $`\zeta <1`$ , equation (59) is similar to the one obtained by Fisher  for the droplet model of condensation, with $`\zeta `$ corresponding to a surface energy exponent. He showed that $`G(x,w)`$ can be expressed as

$$G(x,w)=_0^{\mathrm{}}\frac{x}{x_ce^tx}f(t)𝑑t,$$
(61)
where the auxiliary function $`f(t)`$ vanishes at the origin with an essential singularity of the form $`f(t)\mathrm{exp}(C/t^\zeta ^{})`$, so $`G(x,w)`$ as well as all its derivatives remain finite when $`xx_c`$ from below and $`G`$ has a very weak singularity of the form

$$\mathrm{log}G_{sing}\frac{C}{(x_cx)^\zeta ^{}},$$
(62)
with $`\zeta ^{}=\zeta /(1\zeta )`$. In particular $`\zeta ^{}=1`$ if $`\zeta `$ has the same value $`1/2`$ as for isotropic $`2d`$ percolation.

### 5.2 Numerical results

As discussed in  the phenomenological renormalization approach may still be used to determine the critical fugacity, as the parallel correlation length diverges at the transition. A direct application of the 3-strip PR approach using (35 \- 38) yields results for $`x_c`$ that converge very rapidly to the value expected from (14). As can be seen from table 4 (see the third column), the finite-size estimates $`\stackrel{~}{x}_n`$ of $`x_c`$ are already quite good for narrow strips. To obtain such an accurate numerical check of (14) we had to use a mathematical package allowing for high-precision numbers (the precision of our input data was set up to 30 digits).

The exponent $`\kappa _n`$, however, is found to diverge, signalling a failure of (33) to describe adequately the finite-size scaling of $`\xi _n^{}`$. A more detailed numerical analysis of the behaviour of $`\xi _n^{}`$, for an interaction $`w=8`$ deep in the low-T region, shows that it increases exponentially with the width $`n`$ if one fixes $`x`$ at the critical value $`x_c(w)`$ given by (14):

$$\xi _n^{}(x_c)A\mathrm{exp}(Kn),$$
(63)
and that the fugacity $`\stackrel{~}{x}_n`$ at which $`\xi _n^{}`$ diverges (i.e., $`\lambda _n(\stackrel{~}{x}_n)=1`$) converges exponentially towards $`x_c`$ (see table 4). These results suggest the following finite-size form:

$$\xi _n^{}\frac{C}{(\stackrel{~}{x}_nx)^\lambda }\frac{C}{(x_cx+D\mathrm{exp}(K^{}n))^\lambda },$$
(64)
where $`\lambda =K/K^{}`$. Numerically one finds that $`KK^{}=2.21(2)`$, so we conjecture that $`\lambda =1`$.

The average mass of the clusters

$$<N(x)>=x\frac{\mathrm{log}G}{x}$$
(65)
is dominated by the regular part of $`G(x)`$ and remains finite when $`xx_c`$ from below at fixed temperature, as well as all the higher moments of the mass distribution, so the mass and linear dimensions of the typical clusters also remain finite. This implies that the characteristic transverse length to be compared to the strip width in a scaling expression remains finite at the transition, which is consistent with the observed divergence of the effective exponent $`\kappa _n`$.

In the zero-temperature limit ($`w\mathrm{}`$) the critical fugacity $`x_c0`$, with $`x_cw^21`$, and the generating function $`G(x,w)`$ is dominated by the configurations with the largest possible number of pairs of neighbouring occupied sites at fixed $`N`$. These correspond to compact animals, having no internal holes and the maximum number of independent loops, or equivalently the minimum number of lateral surface sites (i.e., sites with either the smallest abscissa on a given row or the smallest ordinate on a given column of the underlying square lattice). One would expect these animals to be globular rather than elongated, but it is not clear to us how to extract from the transfer matrix results values for the geometrical exponents of the (extremely rare) large clusters in that regime.

Let us note that the present compact animals are different from the ”directed compact percolation clusters” studied by Essam , which just obey the constraint that two diagonally neighbouring sites are necessarily followed by an occupied site. These are themselves a subset of another class of compact clusters studied in the literature, named ”fully directed compact” by physicists  and more precisely ”diagonally convex directed” (DCD) by mathematicians . DCD animals may be characterized as having the minimum accessible perimeter for a given length (i.e., unoccupied sites having an occupied predecessor). They have been exactly enumerated with respect to various parameters  and are expected to be extremely elongated, with $`\nu ^{}=1`$ and $`\nu ^{}=0`$ (where $`\nu ^{}`$ is associated with the local width at a given large distance from the origin, rather than to the overall animal width ).

## 6 Special points

As pointed out in section (2.4), the phase diagram of interacting directed animals in the $`(w,x)`$ plane contains several special points related to various ensembles of animals having specific geometric properties.

* Thus, from (22) the generating function of bond animals is obtained on the line $`w=\frac{2}{(1x)}`$. Using strips of width $`(14,15,16)`$ we find this line intersects the critical line $`x_c(w)`$ for $`w_c=\mathrm{2.285\hspace{0.17em}08}(1)`$. This is in agreement with the value $`y_b=\mathrm{0.285\hspace{0.17em}087\hspace{0.17em}5}(8)`$, obtained for their critical fugacity by Conway et al.  since in the present notations $`w_c=2+y_b`$.
* Bond trees correspond to $`w=2`$ according to (24) and we find that the corresponding $`x_c=\mathrm{0.147\hspace{0.17em}892\hspace{0.17em}24}(1)`$, in agreement with (but more precise than) the critical fugacity $`y_{bt}=2x_c=\mathrm{0.295\hspace{0.17em}785}(10)`$ found in .
* Site trees correspond to the limit $`w0`$ with $`y=wx`$ fixed and finite. They were studied as Model B in  and their critical fugacity was found there to be $`y_{st}=\mathrm{0.368\hspace{0.17em}649}(2)`$.
* Finally for a site fugacity $`x=1`$ one recovers the generating function of strongly embeddable animals, see (27), so their critical fugacity $`y_{emb}=w_c(1)`$. We find $`y_{emb}=\mathrm{0.356\hspace{0.17em}563\hspace{0.17em}26}(1)`$ for this quantity, which does not seem to have been estimated previously in the literature.

A phase diagram summarizing these various results is displayed in figure 2. It is noteworthy that for $`w<w_\theta `$ the plotted quantity $`y_c=wx_c(w)`$, which according to (18) is just the critical fugacity for animals enumerated according to loop number, varies quasi linearly with $`w`$, while this is not at all the case for the analytic continuation of the low-temperature expression (14).

## 7 Crossover region close to the collapse transition

### 7.1 Critical line and crossover exponent

We now focus our attention on the region close to the $`\theta `$-point. As noted above, for $`w>w(\theta )`$ the critical line $`x_c(w)`$ is given by the simple expression (14) and a remarkable consequence  is the absence of a singularity in the bulk free energy per site $`f(T)`$ when $`ww(\theta )^+`$. The collapse transition is just signalled by a singularity in the size-dependent correction term, which may be interpreted as the vanishing of a surface free energy. For $`w<w(\theta )`$ the critical line departs from the analytic continuation of (14), which we denote $`x_D`$ in the following to avoid confusion. The singular part is traditionally written under the form

$$x_cx_DC(w(\theta )w)^{1/\varphi },$$
(66)
where $`\varphi `$ is the crossover exponent.

In most cases it is quite difficult to obtain accurate values for such crossover exponents from numerical data, as generally neither the non-singular part of the critical line nor the value of the critical temperature are known with high precision, making a fit of the form (66) subject to large uncertainties. Here the situation is quite favourable as we can take advantage of the analytic nature of the free energy in the low-$`T`$ phase and of the high-precision knowledge of the critical point $`w(\theta )`$. The remaining difficulty comes from the uncertainties due to the extrapolation procedure used to obtain $`x_c(w)`$, which prevent us from analyzing reliably the data extremely close to $`w(\theta )`$. Fortunately there exists an interval of values of $`w`$ where the extrapolation uncertainties and the corrections to (66) are both small enough to allow a reliable analysis.

The numerical results are displayed as a log-log plot in figure 3 and the slope of the curve yields

$`1/\varphi `$ $`=`$ $`2.53\pm 0.05`$ (67)
$`\varphi `$ $`=`$ $`0.395\pm 0.008.`$ (68)
This corresponds to a weak singularity for the specific heat measured along the critical line, as the associated exponent $`\alpha `$  is negative:

$$\alpha =21/\varphi =0.53\pm 0.05,$$
(69)
in strong contrast with what is found for the $`\theta `$-point of isotropic animals, where $`\alpha 0.48`$ . This also implies that in the vicinity of their collapse transition interacting directed animals should display geometric properties close to those of critical directed percolation clusters up to large scales.

### 7.2 Relation with a percolation exponent

In order to derive a relation between the crossover exponent $`\varphi `$ and a directed percolation exponent, we rewrite (12) under the form

$$\mathrm{log}x_c=\mathrm{log}x_D\underset{N\mathrm{}}{lim}(1/N)\mathrm{log}Pr(N),$$
(70)
where $`x_D=(w2)/w(w1)^2`$. For $`w<w(\theta )`$, the bond probability $`p<p_c`$ and the cluster probability $`Pr(N)`$ is expected to decay exponentially with $`N`$, so the second term on the right-hand side of (70) is no longer negligible when $`N\mathrm{}`$. Equation (70) may then be written

$`(1/N)\mathrm{log}Pr(N)`$ $``$ $`\mathrm{log}(x_c/x_D)\text{for}N1,`$ (71)
$``$ $`{\displaystyle \frac{x_cx_D}{x_D}}\text{for}x_cx_Dx_D,`$ (72)
$``$ $`C[w(\theta )w]^{1/\varphi },`$ (73)
using (66) in the last line. Relation (11) between $`p`$ and $`w`$ still holds in the high-$`T`$ phase, so $`w(\theta )w(p_cp)/(1p_c^2)`$ and we finally obtain

$$(1/N)\mathrm{log}Pr(N)C(p_cp)^{1/\varphi }.$$
(74)
Now, close to $`p_c`$ the probability distribution of cluster sizes is expected to follow a scaling law of the form

$$Pr(N)N^{1\tau }\mathrm{\Psi }[N^\sigma (p_cp)],$$
(75)
for $`N1`$ and $`z=N^\sigma (p_cp)`$ finite, where $`\mathrm{\Psi }(z)`$ is a scaling function finite for $`z=0`$ and we have followed the notations of  for the Fisher exponents $`\tau `$ and $`\sigma `$.

For these two expressions for $`Pr(N)`$ to be compatible it is necessary that

$$\varphi =\sigma .$$
(76)
A very accurate value for $`\sigma `$ may be deduced from Jensen’s recent results  and a standard scaling relation between percolation exponents :

$$\sigma =1/(\beta _{DP}+\gamma _{DP})=\mathrm{0.391\hspace{0.17em}510}(2),$$
(77)
in good agreement with our numerical value (68) for $`\varphi `$, which was obtained without invoking scaling assumptions.

## 8 Conclusion

We have shown that the combined use of the transfer matrix method and of a phenomenological renormalization-group analysis gives very accurate results when applied to interacting directed animals (IDA), in particular close to their collapse transition. The present results are still not quite as precise for the tricritical fugacity and the tricritical exponents as those obtained from the best series calculations for directed percolation , but the method is more flexible so it could also be applied without excessive effort to situations involving more complicated interactions.

Acknowledgments

We thank M. Bousquet-Melou, H. Chaté, B. Derrida and D. Dhar for useful discussions, P. Grassberger for providing several relevant references. One of us (M. K) would like to thank the University P. et M. Curie (Paris 6) for financial support and the members of the Laboratoire de Physique Statistique de l’ENS, where part of this work has been done, for their kind hospitality.