text:
# Domain wall creep in epitaxial ferroelectric Pb(Zr0.2Ti0.8)O3 thin films

## Abstract

Ferroelectric switching and nanoscale domain dynamics were investigated using atomic force microscopy on monocrystalline Pb(Zr0.2Ti0.8)O3 thin films. Measurements of domain size versus writing time reveal a two-step domain growth mechanism, in which initial nucleation is followed by radial domain wall motion perpendicular to the polarization direction. The electric field dependence of the domain wall velocity demonstrates that domain wall motion in ferroelectric thin films is a creep process, with the critical exponent $`\mu `$ close to 1. The dimensionality of the films suggests that disorder is at the origin of the observed creep behavior.

Understanding the propagation of elastic objects driven by an external force in the presence of a pinning potential is a key to the physics of a wide range of systems, either periodic, such as the vortex lattice in type II superconductors Blatter et al. (1994), charge density waves Grüner (1988) and Wigner crystals Andrei and al. (1988), or involving propagating interfaces, such as growth phenomena Kardar (1995), fluid invasion Wilkinsion and Willemsen (1983) or magnetic domain walls Lemerle et al. (1998). In particular, the response to a small external force is of special theoretical and practical interest. It was initially believed that thermal activation above the pinning barriers should lead to a linear response at finite temperature Anderson and Kim (1964). However, it was subsequently realized that a pinning potential, either periodic Blatter et al. (1994) or disordered, Ioffe and Vinokur (1987); Nattermann (1987); Blatter et al. (1994); Chauve et al. (2000), can lead to diverging barriers and thus to a non-linear response, nicknamed creep, where the velocity is of the form $`v\text{exp}(\beta R(f_c/f)^\mu )`$. $`\beta `$ is the inverse temperature, $`R`$ a characteristic energy and $`f_c`$ a critical force. The dynamical exponent $`\mu `$ reflects the nature of the system and of the pinning potential. Despite extensive studies of the creep process in periodic vortex systems Blatter et al. (1994), precise determination of the exponents has proven difficult, given the many scales present in this problem and the range of voltage needed to check the creep law Fuchs et al. (1998). For interfaces, a quantitative check of the creep law has been done recently in ultrathin magnetic films Lemerle et al. (1998), where an exponent $`\mu =0.25`$ has been measured in very good agreement with the expected theoretical value for this system. Quantitative studies of this phenomena in other microscopic systems with other pinning potentials are clearly needed.

In this respect, ferroelectric materials are of special interest. These systems possess two symmetrically equivalent ground states separated by an energy barrier $`U_0`$, as illustrated in Fig. 1 for a tetragonal perovskite structure. Each state is characterized by a stable remanent polarization, reversible under an electric field.

Regions of different polarization are separated by elastic domain walls. The application of an electric field favors one polarization state over the other, by reducing the energy necessary to create a nucleus with a polarization parallel to the field, and thus promotes domain wall motion. In addition to pure theoretical interest, understanding the basic mechanism of domain wall motion in ferroelectrics has practical implications for technological applications, such as high-density memories. In bulk ferroelectrics, switching and domain growth were inferred to occur by stochastic nucleation of new domains at the domain boundary, a behavior observed in BaTiO3 and triglycine sulphate by studies using combined optical and etching techniques Merz (1954); Fatuzzo and Merz (1959); Miller and Weinreich (1960). Early analysis of such motion reported a field dependence of the domain wall speed, $`v\mathrm{exp}[1/E]`$. This bulk system behavior was explained by assuming that domain walls propagate via nucleation along 180 domain boundaries.

In this Letter, we report on studies of ferroelectric domain wall motion in atomically flat single crystal Pb(Zr0.2Ti0.8)O3 films by atomic force microscopy (AFM), allowing noninvasive investigation of domain dynamics with nanometer resolution. In this model system, we identify domain wall motion to be a disorder-controlled creep process. The dynamical exponent $`\mu `$ is found to be close to 1. The activation energy increases significantly from 0.5MV/cm to 1.3MV/cm as the film thickness is reduced from 810Åto 290Å.

The ferroelectric materials investigated were epitaxial c-axis oriented Pb(Zr0.2Ti0.8)O3 thin films, RF-magnetron sputtered onto conducting (100) Nb-doped SrTiO3 substrates Tybell et al. (1998); Triscone et al. (1996). This system allows precise control of film thickness and crystalline quality, exhibiting atomically smooth surfaces with a polarization vector parallel or antiparallel to the c-axis Tybell et al. (1998); Triscone et al. (1996). To study switching dynamics in these films, we used a conductive AFM tip to artificially modify domain structure Tybell et al. (1998); Güthner and Dransfeld (1992); Ahn et al. (1992). Domains were polarized by applying a voltage pulse across the ferroelectric film, between the tip and the substrate. Their sizes were subsequently measured by piezoelectric microscopy Güthner and Dransfeld (1992), as a function of pulse width and amplitude. For every pulse width, we used 12 V pulses to polarize an array of 16 domains, and calculated an average domain size based on their vertical and horizontal radii. The rms error was $``$10%. All domains studied were written in a uniformly polarized area.

Fig. 2a shows the domain radius as a function of pulse width, and three piezoelectric images of ferroelectric domain arrays written with 50$`\mu `$s, 1ms, and 100ms voltage pulses. As can be seen, varying the writing time (the pulse width) markedly changes the size of the AFM-written domains. We observe that domain radius increases logarithmically with increasing writing time for times longer than $``$20$`\mu `$s. Below 20$`\mu `$s, and down to 100ns, the shortest times investigated, domain radius is found to be constant and approximately equal to 20nm, as shown by the shaded area in Fig. 2a Paruch et al. (2001). All of our data suggest that this minimum domain size is related to the typical tip size used for the experiments, whose nominal radius of curvature is $``$20-50nm. In previous work, we have observed that domain size also depends linearly on writing voltage, above a threshold related to the coercive field Paruch et al. (2001). A detailed analysis of the data reveals only well-defined homogeneous domains with regular spacing, as can be seen in the piezoelectric image of a regular 90-domain array, written with 1ms pulses in Fig. 2b. We note that the topographic image of the same area is featureless, with a rms roughness of $``$0.2nm. Within our $``$5nm resolution we do not detect any randomly nucleated domains. The data thus suggest a two-step domain switching process in which nucleation, originating directly under the AFM tip, is followed by radial motion of the domain wall outwards, perpendicular to the direction of polarization.

To analyze this lateral domain wall motion, we note that the force exerted on the wall is given by the electric field $`E`$. To obtain the electric field distribution, we model the tip as a sphere, with a radius $`a`$. The potential on the ferroelectric surface at a distance $`r`$ from the tip is then $`\mathrm{\Phi }\frac{Va}{r}`$, and the local field across the ferroelectric $`E=\frac{Va}{rd}`$, $`V`$ being the applied tip bias and $`d`$ the film thickness. This equation will allow us to relate the change in domain size to the local electric field near the domain boundary. By writing arrays with different pulse widths, and subsequently calculating the average domain size for a given time, we can extract the speed of the domain wall as $`v=\frac{r(t_2)r(t_1)}{t_2t_1}`$ and the corresponding electric field $`E(r)`$ where $`r=(r(t_1)+r(t_2))/2`$. Fig. 3 shows the wall speed as a function of the inverse field for three different film thicknesses. The data fits well to a creep formula

$$v\mathrm{exp}\frac{R}{k_BT}\left(\frac{E_0}{E}\right)^\mu $$
(1)
with $`\mu =1`$. The exact dynamical exponent $`\mu `$ is found to be 1.12, 1.01 and 1.21 for the 290Å, 370Å, and 810Å thick films respectively,with an estimated 10% rms error on the field 111The 10% rms error in the local electric field is given by the precision of the domain size measurements. We used 20nm as the AFM tip radius of curvature.. We find the effective “activation energy” $`\left(R/(k_BT)\right)^{\frac{1}{\mu }}E_0`$ to be 1.321MV/cm, 1.305MV/cm, and 0.506MV/cm for 290Å, 370Å, and 810Å thick films, typically 2 orders of magnitude larger than the applied fields during the polarization process. We note that during AFM writing, the effective field across the ferroelectric is approximately 10 times smaller than the field E 222The exact magnitude of the electric field across the ferroelectric is difficult to quantify, because of a possible gap between the film surface and the tip Hidaka et al. (1997), and the unknown and gradually changing tip shape. Local piezoelectric hysteresis measurements on 120Å- 800Åthick films show that the minimum switching field is $``$16 - 6 times larger than the bulk coercive field. We therefore correct the calculated field E by a factor of 10 to estimate the effective field across the film in our spherical tip model.. This has no effect on the exponent $`\mu `$ governing the exponential velocity dependence, but decreases the effective “activation energy” by a factor of 10, which has been taken into account in the values given above.

Let us now consider the possible microscopic origins of the observed creep behavior. Creep phenomena are a consequence of competition between the elastic energy of a propagating interface, tending to keep it flat, and a pinning potential, preventing it from simply sliding when submitted to an external force. The dynamical exponent in the creep scenario depends on both the dimensionality of the system, and the nature of the pinning potential. Although creep processes are generally associated with the glassy behavior of disordered systems, they can also be observed in a periodic potential if the dimensionality of the object is larger than or equal to 2. For thick films, as for bulk ferroelectrics, the domain wall is a two dimensional object. In this case, from free-energy considerations, and neglecting the anisotropic dipole field present in a ferroelectric, one would expect the exponent to be $`\mu =1`$ Blatter et al. (1994). In PbTiO3, it has been shown theoretically that domain wall energy depends upon whether the wall is centered on a Pb or Ti plane Pöykkö and Chadi (1999); Meyer and Vanderbilt (2001), giving rise to an intrinsic periodic pinning potential. One possible explanation could thus be that the observed creep is due to the motion of the two-dimensional wall in this periodic potential. Note that this scenario is a generalization of the nucleation model developed for bulk ferroelectrics Miller and Weinreich (1960). In order to test this hypothesis we calculated the size of the critical nucleus, using the formula derived by Miller and Weinreich Miller and Weinreich (1960). To estimate $`l^{}`$, the critical length along the c-axis, we used the standard remanent polarization, lattice parameters and dielectric constant values for PZT , the 169mJ/m2 domain wall energy derived for PbTiO3 Meyer and Vanderbilt (2001), and the corrected values for the electric field across the ferroelectric film. We find $`l^{}`$ to vary, depending on the field range, between 200Å and 500Å, 600Å and 1100Å, and 900Å and 1700Å for the 290Å, 370Å and 810Å. The critical nucleus would thus need to be larger than the thickness of the system. Furthermore, the effective “activation energy” calculated for the nucleation model is two orders of magnitude greater than the 0.5-1.3MV/cm determined experimentally. Calculations directly starting from the periodic potentials given in Meyer and Vanderbilt (2001) lead to similar conclusions. These results strongly suggest that the films are in a two dimensional limit, and that the nucleation model or equivalently the motion through a periodic potential does not adequately explain the experimental data.

The creep behavior in the film thus has to be due to disorder in the system, and thereby to the glassy characteristics of a randomly pinned domain wall, with the creep exponent dependent on the nature of the disorder. Defects locally modifying the ferroelectric double well depth $`U_0`$ and giving rise to a spatially varying pinning potential would lead to a “random bond” scenario similar to the one for the ferromagnetic domain walls Lemerle et al. (1998). The exponent $`\mu `$ would be $`\mu =\frac{d2+2\zeta }{2\zeta }`$ where $`\zeta `$ is a characteristic wandering exponent and $`d`$ the dimensionality of the wall. For one dimensional domain walls $`\mu =1/4`$ whereas for two dimensional ones $`\mu 0.50.6`$, hardly compatible with the data. If, however, the defects induce a local field, asymmetrizing the double well, or if there are spatial inhomogeneities in the electric field, the situation is different. For such a “random field” scenario $`\zeta =\frac{4d}{3}`$ Fisher (1986) and thus leads to $`\mu =1`$ for $`1<d<4`$. This scenario is therefore compatible with the observed data. However, further study determining the wandering exponent would be needed to ascertain this precise point 333Only short range elastic forces are considered in the above discussion. Long range interactions would change the elastic energy, and therefore the dynamics of the system, and hence the exponent $`\mu `$. In this case, the random bond scenario could also lead to different exponents and possibly to $`\mu `$ close to 1..

Finally, we note the applications of these results to technological developments. The fact that the domain walls exhibit creep motion with a relatively large exponent $`\mu =1`$ implies a strong stability of ferroelectric domains in thin films, since the induced speed of the domain wall becomes exponentially small as the driving force goes to zero. All domains studied in this work were stable under ambient conditions for the entire 7 day duration of the experiment. Furthermore, sub-$`\mu `$m wide line-shaped domains were stable up to one month Marré et al. (2001). Previously, large, regular arrays of ferroelectric domains with densities of $``$6Gbit/cm2 have been reversibly written Paruch et al. (2001), and densities up to 150Gbit/cm2 have been extrapolated from sizes of individual domains Maruyama et al. (1999). This work identifies the key parameters controlling domain size: the strength, duration, and confinement of the applied electric field. These can be exploited to increase the information storage density in ferroelectric arrays. By using short voltage pulses, for which the resulting domain size is independent of the writing time, and domain-domain separation as small as 10nm, regular arrays with densities of the order $``$30Gbit/cm2 could be written. Fig. 4 shows such an array written on a 370Å thick sample with a density of 28Gbit/cm2. Furthermore, it has been suggested Woo et al. (2001) that the minimum stable domain size is given by the film thickness. Our data, however, show that the minimum stable domain size, $``$40nm, is independent of the thickness used in this study; rather, it is given by the area over which the electric field is applied during the polarization process. Therefore, by confining the field to a smaller area, domain wall creep will be limited, leading to smaller domains and hence even higher information densities.

In conclusion, our studies demonstrate that lateral domain wall motion in ferroelectric thin films is a creep process, governed by a characteristic dynamical exponent $`\mu `$ close to $`1`$. The dimensionality of our films suggests that disorder is at the origin of the observed creep behavior, which inherently explains the measured $`\text{exp}(1/E)`$ dependence of the domain wall speed. The activation energy is found to be around 1MV/cm, decreasing with increasing film thickness. Finally, these results suggest a high degree of stability for ferroelectric domains in low electric fields and identify the key parameters controlling domain size.

###### Acknowledgements.

The authors would like to thank Gianni Blatter, Dima Geshkenbein, and Piero Martinoli for enriching discussions. Special thanks to Steve Brown for carefully reading the manuscript, and to Daniel Chablaix for useful technical developments. This work was supported by the Swiss National Science Foundation through the National Center of Competence in Research “Materials with Novel Electronic Properties-MaNEP” and Division II.


format: text

text:
# Ferroelectric Dynamics in the Perovskite Relaxor PMN

## 1 Introduction - Soft Modes and PNR

Recent neutron scattering studies of the perovskite ($`AB`$O3) relaxors Pb(Mg1/3Nb2/3)O3 (PMN) and Pb(Zn1/3Nb2/3)O3 (PZN) have focused on the low-energy ($`20`$ meV) lattice dynamics to determine whether or not a soft mode picture is relevant to these chemically-disordered systems. pmg\_prl1 ; pmg\_prb ; pmg\_prl2  The answer to this question is not obvious in spite of the similarities between PMN and PZN, and the displacive ferroelectric PbTiO3, which exhibits a phase transition to a tetragonal structure at $`T_c=763`$ K. The occupation of the $`B`$-site by either Mg2+ or Zn2+, and Nb5+ cations creates a local charge imbalance, and gives rise to rapidly varying random fields. ye  Therefore anomalies in the lattice dynamics of both PMN and PZN are to be expected. Of particular interest is the role played by the polar nanoregions (PNR), which for PMN form between 600 and 650 K. burns  Given their polar nature they ought to couple strongly to polar phonon modes, including any soft modes if present, in the relaxors PMN and PZN.

The existence of a soft transverse optic (TO) phonon mode has been documented in numerous perovskites including PbTiO3, BaTiO3, SrTiO3, and KTaO3. In the case of PbTiO3, the frequency of the lowest-lying TO mode drops when cooled from high temperature and condenses at the (first-order) transition temperature $`T_c=763`$ K, transforming the system into a tetragonal state. shirane  Unlike those in other perovskites, the soft mode in PbTiO3 is well defined except for a limited range of reduced wave vectors $`q0.08`$. The zone center $`q=0`$ mode is, however, rather broad and well defined only at temperatures far above $`T_c`$shirane

Pioneering lattice dynamical studies were carried out by Naberezhnov et al. on PMN from 300 K to 900 K. naberezhnov  Fig. 1 shows the TO and TA phonon dispersions for PMN measured at 800 K around (221). A well defined zone center TO mode is observed at this temperature. However, it was identified as a hard TO1 mode in this study because its structure factor was inconsistent with those obtained from neutron diffuse scattering measurements by Vakhrushev et al. vakhrushev95  Prior Raman scattering studies had also found no evidence of a soft mode in PMN. Siny

Subsequent experiments by Gehring et al. were done on PZN and PZN doped with 8% PbTiO3 (PZN-8%PT) for which large good quality single crystals were readily available. Fig. 2 combines a series of constant-$`\stackrel{}{Q}`$ spectra measured on a 4.2 gm single crystal of PZN at 500 K, nearly 100 K above $`T_c`$, into a greyscale contour plot (see reference for experiment details). pmg\_prb  These data show an unusual feature in which the TO branch appears to drop into the TA branch around a $`q`$ of just under 0.15 reciprocal lattice units (1 rlu = 1.545 Å-1), which is a measure of the average size of the PNR. This anomaly, termed “the waterfall” because of its vertical appearance, occurs as a result of a $`q`$-dependent damping in which the TO phonon mode becomes increasingly damped with decreasing wave vector $`q`$, i. e. increasing wavelength. No phonon peaks are observed in constant-$`\stackrel{}{Q}`$ scans for $`q`$ near the zone center, which indicates that these TO modes are overdamped. Thus the waterfall feature is not a true dispersion because the scattering does not come from propagating modes. The waterfall is seen in PMN as well as in solid solutions of both PZN and PMN and PbTiO3 (PZN-$`x`$PT and PMN-$`x`$PT), suggesting it may be common to all relaxor systems. koo ; gop

In comparison to PbTiO3, these results on PZN reveal the presence of a new damping mechanism that effectively freezes out the longest-wavelength TO modes and produces the unusual waterfall feature. The PNR are the most plausible cause of this damping. To test this idea would require heating the PZN crystal to temperatures well in excess of $`T_d760`$ K where the PNR are no longer present. Because this was deemed too close to the decomposition temperature of PZN, this experiment was performed on PMN, for which $`T_d620`$ K as first determined by the optical measurements of Burns and Dacol. burns  As was the case in PZN, an overdamped response was observed in PMN at $`q=0`$ at 500 K. However the data in Fig. 3 show unambiguous evidence of a well defined peak at $`q=0`$ at 1100 K that softens upon cooling. At the same time the mode becomes increasingly damped (broader energy linewidth). Within experimental uncertainty the soft mode becomes overdamped at $`T_d`$, thereby establishing a definite connection between it and the PNR. Equally interesting is the fact that the square of the zone center TO mode energy varies linearly with temperature ($`(\mathrm{}\omega _0^2)(TT_0)`$) as shown in the inset. This is the hallmark of a ferroelectric soft mode. While the softening of this mode is clear, the uncertainties in the data obtained from this small 0.1 cm3 sample do not permit a reliable determination of $`T_0`$, which should be close to 400 K. viehland  Data obtained on a much larger crystal (0.4 cm3) are discussed in a later section and are combined with the data shown in Fig. 3 (see Fig. 5).

## 2 The Connection to the PNR

If the PNR are in fact the result of the condensation of the soft TO mode in PMN, then the ionic displacements that characterize the diffuse scattering intensities below $`T_d`$ must be consistent with those that characterize the soft mode above $`T_d`$. In particular, since the ionic motions associated with the zone center TO vibrational mode cannot shift the center of mass of the PMN unit cell, the static displacements observed below $`T_d`$ must also satisfy this constraint. But they do not. Vakhrushev et al. measured the neutron diffuse scattering intensities near 16 reciprocal lattice points on a single crystal of PMN, and obtained the following relative ionic displacments $`\delta `$ (normalized such that $`\delta _{\mathrm{Pb}}=1`$): $`\delta _{\mathrm{MN}}=0.18`$, $`\delta _{\mathrm{O}(1)}=0.738`$, $`\delta _{\mathrm{O}(2)}=0.549`$vakhrushev95  These displacements do not satisfy the center-of-mass constraint. Yet if one believes that the static and dynamic ionic displacements are related to the same polar state, then one is faced with an unavoidable discrepancy that must be resolved.

This apparent contradiction was recently resolved by Hirota et al. who noticed that if a uniform shift of the proper magnitude were subtracted from each of the ionic displacments, then the remaining displacements satisfied the center-of-mass constraint. In other words, for each ion they assumed $`\delta _i=\delta _i^{\mathrm{cm}}+\delta _{\mathrm{shift}}`$, where $`i`$ is the ionic index. Using this model, the size of the uniform shift was determined to be $`\delta _{\mathrm{shift}}=0.58`$, assuming an average oxygen displacement of $`\delta _\mathrm{O}=0.64`$. The corresponding values for the $`\delta _i^{\mathrm{cm}}`$ are: $`\delta _{\mathrm{Pb}}^{\mathrm{cm}}=0.42`$, $`\delta _{\mathrm{MN}}^{\mathrm{cm}}=0.40`$, $`\delta _\mathrm{O}^{\mathrm{cm}}=1.22`$hirota

Calculations of both the inelastic and diffuse scattering structure factors were made using these values. Fig. 4 (a) shows the calculated intensities of the TO phonon mode at the (101), (200), and (300) reciprocal lattice points as a function of the ratio $`S=S_2/S_1`$ of the Slater ($`S_1`$) and Last ($`S_2`$) modes. harada  Inelastic scans measured at (200) and (300) (see Fig. 6) fix the value of $`S`$ close to 1.5, consistent with that for PbTiO3. Fig. 4 (b) shows the corresponding diffuse scattering intensities for this value of $`S`$ as a function of the single parameter $`\delta _{\mathrm{shift}}`$. Agreement with experiment, which indicates extremely little diffuse scattering around (200), is obtained precisely at $`\delta _{\mathrm{shift}}=0.58`$. Because of the uniform shift $`\delta _{\mathrm{shift}}`$, the PNR are said to result from a “phase-shifted condensed soft mode.” The shift of the PNR is along their respective polar directions.

## 3 Ferroelectric Dynamics in PMN below $`T_c`$

Above $`T_d620`$ K PMN exhibits dynamical properties characteristic of a normal ferroelectric in its cubic phase that are surprisingly similar to those observed in PbTiO3. The formation of PNR below $`T_d`$ appear to stifle the propagation of long-wavelength polar TO modes due to their finite size. The concurrent development of diffuse scattering intensity at $`T_d`$, first observed by Naberezhnov et al. naberezhnov , is consistent with the picture of PNR condensing from the soft TO mode. Yet PMN is reported to remain cubic on average at all temperatures (at least down to 5 K). mathan  On the other hand, it is known that a macroscopic ferroelectric phase can be induced in PMN by cooling it in a small electric field. After removal of this field, the ferroelectric state in PMN is lost via a first-order transition when heated above $`T_c=213`$ K. ye2 ; ye3  Therefore the lattice dynamical studies of PMN were extended to low temperature, placing special emphasis on the temperature region near $`T_c`$.

Fig. 5 summarizes recent neutron inelastic scattering results obtained by Wakimoto et al. on a large 0.4 cm3 single crystal between 10 and 750 K. wakimoto  By far the most remarkable feature is the recovery and subsequent behavior of the soft mode where $`(\mathrm{}\omega _0)^2(T_0T)`$ below $`T_c`$ shown in Fig. 5 (b) that is typical of many ferroelectrics in their ordered states. This implies that a phase transition from the relaxor state to a low-temperature, short-range ordered ferroelectric state, occurs at $`T_c`$ in the absence of an applied field. But this contradicts the average cubic symmetry of PMN. However, if the length scale of the ferroelectric order were $`100`$ Å, then one could plausibly resolve this contradiction since this length scale would appear as long-range order in neutron phonon measurements, but short-range order to x-rays. Indeed, regions of ferroelectric order $`100`$ Å in size have been reported by de Mathan et almathan  Hypersonic damping Tu  and Raman scattering Siny  studies have also reported anomalies at or near $`T_c`$, as has a change in the topology of the PNR, Vakhrushev , all of which lend support to the idea of a low-temperature ferroelectric state in PMN. It is interesting to note that in a zone center mode is observed at low (25 K) temperature in PZN at essentially the same energy as that found in PMN (10.5 meV). pmg\_prb

Fig. 5 (a) shows an anomalous damping of the TA mode that coincides with the overdamping of the soft mode and the appearance of the PNR at $`T_d`$. This suggests a non-uniform distortion of the lattice which then vanishes below $`T_c`$, ostensibly when the short-range ordered ferroelectric phase is established. A similar TA broadening has also been observed in PMN-20%PT by Koo et al. koo

## 4 Mode Coupling Analysis

The scattering cross section used to fit these inelastic data is based on the coupled-mode description first used by Harada et al. to fit asymmetric phonon lineshapes observed in BaTiO3harada2  Typical data measured in two different zones, (200) and (300), and taken above $`T_d`$ at 690 K, are shown for PMN in Fig. 6. The coupled-mode cross section is parametrized by a coupling-constant $`\lambda `$, along with inelastic structure factors $`F_1`$ and $`F_2`$ for the coupled TO and TA modes, respectively. Additional parameters are the phonon energies $`\mathrm{}\omega _1`$ and $`\mathrm{}\omega _2`$, and the linewidths $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$. All of these parameters were allowed to float freely between the two scans. The values for the TA and TO phonon energies were identical to within 0.1 meV between (200) and (300), as required, with $`\mathrm{}\omega _1=10.8`$ meV, $`\mathrm{}\omega _2=4.2`$ meV, and represent a useful cross check of the integrity of the data. However the ratio of the dynamic structure factors $`F_1/F_2`$ is markedly different in each zone. At (200) the ratio is 2.3, whereas at (300) it is 27.5, nearly an order of magnitude different. The coupling constant $`\lambda `$ has nearly the same absolute value in each zone (-14.1 versus 17.6 for (200) and (300) respectively), but opposite sign. Physically this is due to a sign change in the TO structure factor $`F_1`$, not $`\lambda `$.

Fig. 7 conveys the manner in which the coupled-mode cross section varies with the TO mode linewidth $`\mathrm{\Gamma }_1`$ (corresponding to an increased damping), as well as with a change in $`\lambda `$. This cross section gives an excellent model description, shown in panel (a), of the $`q`$-dependent damping observed in both PZN and PMN below $`T_d`$. Two cases are presented in panels (b) and (c) following Bullock et al. bullock  Panel (b) corresponds to well-separated TA and TO modes with energies of 5 and 10 meV, respectively, and identical linewidths of 1 meV. Panel (c) corresponds to TA and TO modes with energies of 6 and 8 meV, and substantially different linewidths of 0.5 and 4 meV, respectively, such that there is substantial spectral overlap. bullock  In both cases there is a clear transfer of spectral weight from the TO mode to the TA mode when the mode coupling is turned on ($`\lambda 0`$).

## 5 Conclusions

We have given a brief overview of the lattice dynamics of PMN and PZN. In PMN a well-defined soft mode is observed at 1100 K that becomes increasingly damped upon cooling, and eventually becomes overdamped at $`T_d`$. The overdamping of the soft mode and the appearance of the diffuse scattering both at $`T_d`$ strongly point to a direct connection between the soft mode and the PNR. The resolution of the discrepancy between the static and dynamic structure factors using the “phase-shifted condensed soft mode” model of Hirota et al. provides a convincing framework in which this connection can be understood. The damping of the soft mode that is present at temperatures above $`T_d`$ is believed to be a result of the fact that the PNR exist as dynamical entities at higher temperature.

The soft mode dynamics of PMN below $`T_c`$ are characteristic of an ordered ferroelectric phase, and stand in contrast to the average cubic structure reported for PMN. The assumption of short-range order ($`100`$ Å) can be used to reconcile these two sets of experimental observations. At all temperatures, a coupled-mode cross section, properly convolved with the instrumental resolution function, provides an excellent description of the neutron inelastic data.

A very different interpretation of the soft mode dynamics in PMN has recently been proposed by Vakhrushev and Shapiro vs  in which the existence of a soft “quasi-optic” (QO) mode, distinct from the zone center TO mode, is claimed to explain the observed Curie-Weiss behavior of the dielectric susceptibility above $`T_d`$. The QO-mode is derived from a mode-coupling analysis of neutron inelastic lineshapes measured between 490 and 880 K, and is attributed to the intrinsic chemical disorder of PMN. More studies of PMN at higher resolution are clearly needed to reconcile these different points of view.

We thank A. A. Bokov, K. Hirota, V. Kiryukhin, T. -Y. Koo, K. Ohwada, S. M. Shapiro, C. Stock, N. Takesue, S. B. Vakhrushev, and H. You for stimulating discussions. Financial support from the U. S. DOE under contract No. DE-AC02-98CH10886, and the Office of Naval Research under Grant No. N00014-99-1-0738, is acknowledged. Work at the University of Toronto is part of the Canadian Institute for Advanced Research and is supported by the Natural Science and Engineering Research Council of Canada. We acknowledge the support of the NIST Center for Neutron Research, the U. S. Department of Commerce, for providing some of the neutron facilities used in the present work.


text:
# Are high-redshift damped Lyman𝛼 galaxies Lyman-break galaxies? 1footnote 11footnote 1Based on observations made with the NASA/ESA Hubble Space Telescope. 2footnote 22footnote 2Based on observations collected at the European Southern Observatory, Paranal, Chile (ESO Programme 63.O-0618)

## 1 Introduction

At any redshift the majority of the neutral hydrogen in the Universe is contained in the absorbers of highest column density, the damped Ly$`\alpha `$ (DLA) absorbers. Observations of the column density distribution functions $`f(N,z)`$ for both HI and metals can be used to compute the cosmic mean densities $`\mathrm{\Omega }_{\mathrm{HI}}`$ and $`\mathrm{\Omega }_m`$ as functions of redshift $`z`$. The metallicities of the DLA systems have a large scatter at all redshifts (Pettini et al. 1999; Prochaska & Wolfe 1999); the mean metallicity is 5-10% solar at $`z23`$ and, within the current statistical uncertainties, may or may not increase with decreasing redshift, in accordance with models of cosmic chemical evolution (Kulkarni & Fall 2002). Selection biases caused by dust could dilute the apparent rate of chemical enrichment in the DLA absorbers (Pei & Fall 1995; Boissé et al. 1998). The dust-to-metals ratio also has a large scatter among individual DLA absorbers, but the mean value appears to remain roughly constant with redshift at about the value in the Milky Way (Pei, Fall, & Bechtold 1991; Pettini et al. 1997). The presence of metals implies that DLA absorbers are, or have been, in a state of active star formation and it is therefore interesting to ask whether those stellar populations are themselves visible, and if they are related to other known classes of high redshift objects. In this paper we shall narrowly focus on a single question: “Are high-redshift damped Ly$`\alpha `$ galaxies Lyman-break galaxies?”.

The Lyman-break galaxies (LBGs, Steidel et al. 1996) are starburst galaxies at high redshift, recognisable by their relatively flat restframe ultraviolet continuum, and a sharp discontinuity at the Lyman limit. Under the definition of LBGs we include all objects identified by their characteristic continuum shape as such, on the basis of broadband photometry alone i.e. without requiring spectroscopic confirmation. The typical brightness of an LBG is not a well defined quantity as it depends strongly on the flux limit of the survey. In particular one may consider two distinct types of LBG survey, where surveys concerned solely with the determination of photometric redshifts have significantly fainter flux limits than those which classify as LBGs only galaxies with spectroscopically confirmed redshifts. It is of interest to note that to an AB magnitude limit $`H=26.5`$, i.e. flux limited in the restframe optical continuum, nearly all galaxies in the Hubble Deep Field N in the redshift range $`2<z<3.5`$ are classified as LBGs (Papovich, Dickinson, & Ferguson 2001). With the caveat that this sample is small, it may be concluded that for observations in the optical and near-infrared, to very faint magnitudes, any galaxies with steep continuum, either due to reddening or to the presence of a population of old stars, are a minority population. So at these wavelengths, and down to the AB magnitude limit $`H=26.5`$, the term ‘Lyman-break galaxy’ is essentially synonymous with ‘high-redshift galaxy’.

The global rate of star formation at these redshifts measured from emission from the LBGs is consistent with the rate of star formation inferred from the change with cosmic time of $`\mathrm{\Omega }_{\mathrm{HI}}`$ and $`\mathrm{\Omega }_m`$ in the DLA absorbers (Pei, Fall, & Hauser 1999), modulo uncertain corrections in both calculations for the effects of dust. The simplest interpretation of these results is that the DLA absorbers are the reservoirs of gas from which the stars in LBGs are forming. If this is the case, for every high-redshift DLA absorber detected in the spectrum of a quasar there should be stellar emission visible from a Lyman-break galaxy at the absorption redshift, coincident with or near to the line of sight to the quasar 333The covering factor of DLA absorbers is greater than the fraction of sky covered by detectable optical emission from LBGs. So the absorbing gas would extend beyond the optically visible stellar emission.. Yet searches for this stellar emission have had very little success (for reviews and discussions of past surveys see Møller & Warren 1993; Kulkarni et al. 2000). The large number of unsuccessful attempts to identify the galaxy counterparts of DLA absorbers, hereafter “DLA galaxies”, indicates that most DLA galaxies at $`z>2`$ have very small impact parameters or are very faint (e.g. Steidel, Pettini, & Hamilton 1995a). Indeed searches that have produced lists of candidates pending spectroscopic confirmation (e.g. Ellison et al. 2001a), as well as the very limited set of confirmed identifications available at present (Fynbo, Møller, & Warren 1999), both suggest that DLA galaxies mostly are fainter than $`L^{}`$ LBGs, a result supported by upper limits on H$`\alpha `$ emission (Bunker et al. 1999; Kulkarni et al. 2001). This is seen by some to contradict the view that DLA galaxies are Lyman-break galaxies. In particular the theoretical prediction that the angular momentum distribution of dark matter halos should cause LBGs and DLA galaxies to form two “quite distinct populations” (Mo, Mao, & White 1999) has received much attention, and it has been suggested that DLA galaxies make up a separate low surface brightness galaxy population (Jimenez, Bowen, & Matteucci 1999).

At redshifts $`z0`$ the DLA galaxies are much easier to identify, but it is far from clear how they are related to high-redshift DLA galaxies. Nevertheless it has been found that also the low redshift DLA galaxies typically are sub $`L^{}`$ (Steidel et al. 1995b; Lanzetta et al. 1997; Miller, Knezek, & Bregman 1999; Cohen 2001; Turnshek et al. 2001; Bouché et al. 2001).

In this paper we present results from deep images of three spectroscopically confirmed DLA galaxies observed with the HST instruments STIS and NICMOS, and make a comparison of their emission properties with the emission properties of LBGs. The images presented here are the first combined (both STIS and NICMOS) results from a large campaign of HST observations searching for 18 DLA galaxies towards 16 quasars. The targets, the strategy, and the goals are presented in detail in our first paper (Warren et al. 2001) which also presented the results of the NICMOS imaging campaign. Follow-up spectroscopy is being carried out with the ISAAC and FORS instruments on the ESO VLT. In §2 we present a discussion of how to use such a search for DLA galaxies to test the hypothesis that DLA galaxies are LBGs. In §3 we describe briefly the HST observations, and the measurements of the images. In §4 we make a comparison of the measured properties, including sizes, colours, and morphologies, with the properties of LBGs of similar absolute magnitude. In §5 we address the question posed at the top, and present our conclusions. Cosmological parameters of $`\mathrm{\Omega }_{}=0.3`$, $`\mathrm{\Lambda }_{}=0.7`$, $`h=\mathrm{H}_{}/100`$ are assumed, and the term high redshift implies in this paper $`2<z<3.5`$.

## 2 Interpreting searches for DLA galaxies

Because of the glare from the quasar the detectability of a DLA galaxy is dependent both on the brightness of the galaxy and the impact parameter (in addition to the brightness of the quasar and the stability of the point spread function (psf)). Any sample of DLA galaxies, including the small sample of three discussed here, is biased towards brighter objects and larger impact parameters. Yet if DLA galaxies are LBGs, the measured emission properties of any DLA galaxy must fall within the range of emission properties for LBGs of the same absolute magnitude. This statement is true regardless of how the DLA galaxies were selected. This is the basic prediction that we test in this paper. We compare the measured half-light radii, and colours of DLA galaxies, as well as the morphologies, and Ly$`\alpha `$ emission properties, against the same properties measured for LBGs of similar absolute magnitude.

If DLA galaxies are LBGs they will nevertheless have a different luminosity distribution compared to a flux limited sample of LBGs, because DLA absorbers are detected in proportion to their gas cross section. It follows that the luminosity distribution of DLA galaxies is given by the LBG luminosity function weighted by the luminosity dependence of the gas cross section. A simple calculation illustrates the importance of this for the interpretation of searches for DLA galaxies. For a luminosity function of Schechter form, and a power-law relation between gas radius and luminosity $`RL^t`$, the mean luminosity of galaxies selected by gas cross section is $`\overline{L}=L^{}\mathrm{\Gamma }(2+\alpha +2t)/\mathrm{\Gamma }(1+\alpha +2t)`$ (Fynbo, Møller, & Warren 1999; see also Wolfe et al. 1986; Impey & Bothun 1989). For LBGs at $`z3`$, Adelberger & Steidel (2000) measure $`\alpha =1.57`$. Depending on the value of $`t`$, the value of the mean luminosity can be much less than $`L^{}`$. For example for $`t=0.4`$, $`\overline{L}=0.23L^{}`$. In fact there is some theoretical support for a value of $`t0.4`$ (Haehnelt, Steinmetz, & Rauch 2000), based on the kinematics of DLA absorbers as measured from the absorption profiles of low ionisation species. So the notion that DLA galaxies are typically fainter than $`L^{}`$ is perfectly compatible with the hypothesis that DLA galaxies are Lyman-break galaxies. We note that a consequence of this is that the value of any property which depends on galaxy luminosity will differ between galaxy samples selected either by absorption (DLA galaxies) or by emission (LBG galaxies). For example, since spatial clustering is weaker for less massive (hence less luminous) galaxies (the effects of natural biasing in a hierarchical clustering scenario), we expect DLA galaxies to have a lower clustering amplitude than the emission-selected LBGs.

In the longer term, when we have completed our HST survey, we will have a sample of DLA galaxies complete within well defined selection criteria. Using the measured impact parameters we will be able to determine the form of the relation between radius and luminosity. In fact we have already attempted this (Fynbo et al. 1999). One point to note is that the conclusion that DLA galaxies are much fainter than $`L^{}`$ is quite robust, and insensitive to the uncertainty of the faint-end slope $`\alpha `$ of the Lyman-break galaxy luminosity function. For a different value of $`\alpha `$, the value of $`t`$ must be adjusted in order to reproduce the observed small impact parameters of DLA galaxies. The result is that the value of $`\overline{L}`$ changes little, because the exponent $`\alpha +2t`$ appearing in the expression for $`\overline{L}`$ remains approximately constant.

## 3 Properties of three DLA galaxies

In this section we detail the measurement of the emission properties of the three DLA galaxies, which will be compared against the LBGs in the next section.

### 3.1 Observations and targets

The larger imaging campaign is described by Warren et al. (2001). We are obtaining high-resolution optical and near-infrared images of quasars, searching for faint galaxies close to the quasar line of sight. These are candidate DLA galaxies, targets for spectroscopic follow up. In Warren et al. (2001) we list the 16 target quasars and their coordinates, we tabulate details of the 18 DLA absorbers, present details of the NICMOS observations, and the list of NICMOS candidates. Each quasar was observed for three orbits with the NIC2 camera and the F160W filter. The NICMOS images reach to an AB magnitude typically of $`H_{160}=25`$. A total of 41 candidate counterparts were detected in boxes of side $`7.5\mathrm{}`$ centered on each quasar. Spectroscopy of this candidate list is $`25\%`$ complete. So far one new spectroscopically confirmed DLA galaxy has been discovered. The other two DLA galaxies discussed here are previous discoveries.

In Møller et al. (in preparation) we will present the candidates detected using STIS. With STIS we are imaging in the 50CCD configuration i.e. without a filter in order to reach as deep as possible. In this configuration the effective central wavelength is 5851.5Å and the FWHM is 4410.3Å so 50CCD corresponds to a very wide V magnitude. Each quasar was observed for two orbits, at different orientations, resulting in a detection limit on the combined image of typically $`V_{50}=27`$ (AB magnitude). Because the STIS images reach much fainter magnitudes than the NICMOS images, the number of candidates in the STIS images is typically a factor 2-3 higher than in the NICMOS images.

Fig. 1 reproduces the STIS images of the fields of the three quasars PKS 0528–250, Q 2206–1958, and 2233.9+1318, where the frames from both orbits have been registered and summed. The images show the field after subtraction of the quasar image. Strong residuals near the center of the quasar image have been masked. The DLA galaxies are indicated in each frame, labeled using the numbering scheme in the NICMOS paper. The measured impact parameters are $`1.14\mathrm{}`$, $`0.99\mathrm{}`$, $`2.51\mathrm{}`$, respectively. In each case the confirmed galaxies are the candidates in the NICMOS frames closest to the line of sight to the quasar. The three DLA galaxies have all been confirmed by spectroscopic detection of Ly$`\alpha `$ in emission, but were originally found with three different search techniques. N-7-1C was identified on a deep Ly$`\alpha `$ narrow band image obtained from the ground, N-14-1C was found on the image from our current STIS campaign (Fig. 1), and is also detected in the NICMOS image, and N-16-1D was found on ground based images with the Lyman-break technique. Careful psf-subtraction was in all cases instrumental in the discovery. We now briefly summarise previous results obtained on each of the three DLA galaxies.

#### 3.1.1 N-7-1C

The galaxy N-7-1C was discovered by Møller & Warren (1993), who called it S1. The DLA absorber lies at a similar redshift to the quasar. Using spectroscopic data Warren & Møller (1996) and Møller, Warren, & Fynbo (1998) argued that the physical separation is sufficiently large that the ionising flux from the quasar is not important. This result was supported by Ge et al. (1997) who concluded that the distance between the DLA absorber and the quasar must be larger than 1 Mpc, on the basis of ionization modelling (see also Ledoux et al. 1998; Ellison et al. 2001b). Møller & Warren (1998) obtained WFPC2 images of N-7-1C and showed that the luminosity profile is similar to the profile of LBGs at similar redshift. Ge et al. (1997) reported a metallicity of the DLA absorber of 10% solar and a dust-to-gas ratio of 8% of the Milky Way value. Lu, Sargent, & Barlow (1997) were not able to fit the absorption spectrum of this DLA absorber into the rotating-disk model of Prochaska & Wolfe (1997, 1998), which requires a ‘leading-edge asymmetry’.

#### 3.1.2 N-14-1C

The galaxy N-14-1C is a new discovery. The spectrum confirming the identification is provided in Fig. 2. Details are presented in Tables 1 and 2. The metallicity of the corresponding absorber is 1/3 of the solar value, but despite being one of the most metal rich DLA absorbers known there is no evidence for dust (Prochaska & Wolfe 1997). In this case Prochaska & Wolfe (1997) found evidence for the leading edge asymmetry in the low ion absorption lines, which they interpret as evidence for a rotating disk.

#### 3.1.3 N-16-1D

The galaxy N-16-1D was first reported by Steidel, Pettini, & Hamilton (1995a) who called it N1. Spectroscopic confirmation of the redshift was obtained by Djorgovski et al. (1996). Lu, Sargent, & Barlow (1997) reported a metallicity ‘typical of DLA galaxies at such redshifts’ (\[Fe/H\] = -1.4). They also found evidence for leading-edge asymmetry in the absorption spectrum.

### 3.2 Results

Photometry and profile fitting. In Table 1 we summarise the results of STIS and NICMOS photometry, and profile fitting of the galaxies. The columns list successively the quasar name, the galaxy name, the STIS and NICMOS AB magnitudes, the galaxy impact parameter $`b`$, i.e. the angular offset of the galaxy from the quasar, and the details of the deconvolved surface-brightness profiles measured on the STIS frames. The measurement of N-14-1C was problematic, and is described later. For N-7-1C and N-16-1D we followed the same procedure described by Warren et al. (2001). Briefly, the best fit galaxy profile is determined by convolving parameterised models with the point spread function, and searching for the $`\chi ^2`$ minimum. The $`68\%`$ confidence intervals on each parameter, marginalising over the other parameters, are computed on the basis of the Poisson errors for each pixel by locating the orthogonal tangent planes to the $`\mathrm{\Delta }\chi ^2=1.0`$ ellipsoid. We used the Sersic model, where the surface brightness as a function of radius $`r`$ is $`\mathrm{\Sigma }=\mathrm{\Sigma }_{1/2}\mathrm{exp}\{B(n)[(r/r_{1/2})^{1/n}1]\}`$, where $`r_{1/2}`$ is the deconvolved half-light radius. Note that $`n=1`$ corresponds to an exponential profile, and $`n=4`$ to a de Vaucouleurs profile (for details see Warren et al. 2001). Table 1 lists $`r_{1/2}`$, the value of $`n`$, the galaxy ellipticity, and the parameter $`\varphi `$ which is the angle between the galaxy major axis and the line pointing to the quasar. For these two sources the STIS magnitudes are total magnitudes, calculated by extrapolating the model fits to infinite radius. Because of the low $`S/N`$ of the NICMOS frames total NICMOS magnitudes were measured by determining the scaling of the STIS models that provided the best fit to the NICMOS data.

The galaxy N-14-1C has an irregular morphology (Fig. 1) with a bright knot at the end of a diffuse structure. The measured parameters of the fit varied greatly with the size of the fitting region. Therefore we resorted to aperture photometry to measure magnitudes, with a circular aperture of 0.6″ diameter. In order to obtain an estimate of the half-light radius we employed a model with fewer parameters, a circularly symmetric exponential model, which was much less sensitive to box size. Since the model is not a good fit we assigned a large uncertainty, $`0.1\mathrm{}`$, to this measurement.

Spectroscopy. As stated above the galaxy N-14-1C is a new discovery. The confirming spectrum was obtained with the FORS1 instrument on the European Southern Observatory 8.2m telescope UT1 on the night of 1999 August 12. The integration time was 8000s. With a slit width of $`1.3\mathrm{}`$ and the 600B grating the resolution was 3.0Å. We used the SPSF method of Møller (2000) to extract the galaxy spectrum. The software optimises the spectral extraction of a faint galaxy close to a bright companion (here the quasar). In Fig. 2 we show the region of the spectrum around the wavelength of the DLA absorber, $`z_{DLA}=1.9205`$. An emission line at a wavelength of 3553.4Å, corresponding to Ly$`\alpha `$ at redshift of $`z=1.9229\pm 0.0006`$, is clearly visible. The close agreement between the emission and absorption redshifts, and the small impact parameter of the galaxy confirm that N-14-1C is the DLA absorber galaxy counterpart.

In Table 2 we summarise relevant properties of the Ly$`\alpha `$ emission line of the three galaxies. In successive columns are listed the galaxy name, the quasar redshift, the column density and redshift of the absorption line, and then the following properties of the emission lines: redshift, velocity difference relative to the absorber, FWHM, flux, and rest-frame equivalent width.

## 4 Comparison with other high-redshift galaxies

In this section we compare the measured emission properties for the three DLA galaxies against the distribution of emission properties measured for LBGs, also observed with HST. We begin by considering the half light radii and the observed optical-to-near-infrared colours. To obtain a comparison sample of LBGs we started with the catalogue of Cohen et al. (2000) of redshifts for galaxies in the Hubble Deep Field (HDF) N. We selected the 26 galaxies with spectroscopic redshifts in the range $`2.0<z<3.5`$. This sample was reduced to 24 by removing two galaxies where the published redshifts are the subject of debate (Papovich et al. 2001). Marleau & Simard (1998) have measured the surface brightness profiles of galaxies in the HDF using Simard’s GIM2D software. We matched 21 of the 24 galaxies with their catalogue, and this forms our comparison sample of LBGs.

Half-light radii: In the list of Marleau & Simard (1998) 16 of the 21 galaxies show a single component. The other five galaxies are decomposed into a total of 14 components. We extracted from their catalogue the half light radius $`r_{1/2}`$ and the total apparent AB magnitude for the WFCAM F606W passband, $`V_{606}`$, for each of the 30 galaxies or components. (Note that the pivot wavelength for the WFCAM F606W configuration is very close to the pivot wavelength for the STIS 50CCD passband.) These data are presented in Fig. 3 (upper plot), plotting half-light radius against absolute magnitude, $`M_{1500}`$ (the AB absolute magnitude at a rest wavelength $`\lambda =1500`$Å – we computed the small $`k`$ corrections assuming a continuum shape $`f_\nu \nu ^{0.5}`$). The LBGs are plotted as open symbols. Our measurements of the same parameters for the three DLA galaxies are plotted as filled symbols. There is no strong trend of size with luminosity for the LBGs. The DLA galaxies have similar luminosities to the faintest of the LBGs plotted. The range of half-light radii for the DLA galaxies is similar to the range for the LBGs, although the half-light radius for N-14-1C is a little larger than the LBG radii plotted. Note, however, that it is quite likely that Marleau and Simard would have decomposed a galaxy like N-14-1C into two components.

Marleau & Simard (1998) found that the LBGs mostly have small bulge-to-disk ratios. The small values of the Sersic parameter $`n`$ listed in Table 1 for N-7-1C and N-16-1D imply similarly disk-dominated profiles. At present we cannot quantify how N-14-1C fits into this pattern.

Optical-to-near-infrared colours: Papovich et al. (2001) provide optical-to-near-infrared ($`V_{606}H_{160}`$) colours for high-redshift galaxies in the HDF N. For the same 21 galaxies considered above they list 23 components. Since there is not a one-to-one correspondence between their list of components and the list of Marleau & Simard (above) we have used the total apparent magnitudes from Papovich et al. to compute absolute magnitudes (their total magnitudes are typically 0.3 mag. fainter than the corresponding values in Marleau & Simard 1998).

In Fig. 3 (lower plot) we plot $`VH_{160}`$ colour against absolute magnitude for the LBGs and DLA galaxies, using the same symbols as in the upper plot. As noted by Papovich et al. (2001) there is a marked trend of colour with luminosity, in the sense that galaxies of lower luminosity are bluer. They point out that this trend could reflect a sequence in age, in metallicity, or in dust content, but they do not draw any conclusions on this point. The DLA galaxies are very faint and blue, and have colours in the same range as LBGs of the same absolute magnitude.

Morphology: In addition to the size/luminosity relation, one may consider the more detailed morphology of the objects. Are DLA galaxies for example more disturbed, or do they have more components than LBGs? To address this question we searched the list of Fontana et al. (2000) for LBGs in the HDF S with redshift and V magnitudes as similar as possible to each of the three DLA galaxies. To be specific we selected, for each DLA galaxy, the four LBGs with $`V_{606}`$ and $`z_{phot}`$ that gave the smallest value of the sum $`(V_{606}V_{50,DLA})^2+(z_{phot}z_{DLA})^2`$. We then extracted a $`2\mathrm{}\times 2\mathrm{}`$ image centred on each object. The four objects are shown in Fig. 1 next to the corresponding DLA galaxy. The three DLA galaxies themselves span a wide range in general morphology. N-7-1C has a regular single core morphology, N-14-1C is dominated by a single bright knot at the end of an extended structure, and N-16-1D is elongated with two or three compact components. Qualitatively, one can find examples that resemble each of the DLA galaxies amongst the Lyman-break galaxy images.

Ly$`\alpha `$ emission/absorption: The restframe line equivalent widths (EWs) for the Ly$`\alpha `$ emission lines from the DLA galaxies N-7-1C and N-16-1D are given in the literature as 51 and 37 Å, but with large errors due to uncertainties in determining their broad band fluxes. We have here calculated EWs for all three DLA galaxies as follows: The Ly$`\alpha `$ line falls in all cases inside the 50CCD passband, so we first subtract the measured Ly$`\alpha `$ flux from the $`V_{50}`$ flux. A linear interpolation of the $`H_{160}`$ and $`V_{50}`$ (corrected for Ly$`\alpha `$ ) fluxes then allow a computation of the continuum flux at the wavelength of the redshifted Ly$`\alpha `$ line. The resulting EWs are listed in Table 2, and cover the range 63 to 83Å. The median EW found here is larger than that for LBGs (about 50% of the LBGs have Ly$`\alpha `$ in absorption so their median EW is close to zero). The distribution of Ly$`\alpha `$ EWs for confirmed DLA galaxies will inevitably be skewed towards large emission line EWs as this makes spectroscopic confirmation easier. Nevertheless the values reported here lie within the measured range for LBGs (EWabs = 70Å to EWem = 100Å in a field 92% completed for objects having $`R25.0`$, Steidel et al. 2000), and so are in agreement with the basic prediction we are testing.

Table 2 also lists the velocity difference between the Ly$`\alpha `$ emission and absorption lines $`\mathrm{\Delta }v`$ (Ly$`\alpha `$ –DLA) for the three DLA galaxies. Part of the velocity difference could be due to relative motion, either in the form of random motion of sub-clumps, of ongoing merging processes (Haehnelt, Steinmetz, & Rauch 1998) or of rotation (Prochaska & Wolfe 1997). However, it is noticeable that in all cases the emission line is redshifted relative to the DLA absorber. The lines are also rather broad $`300600`$ km s-1 FWHM. The same characteristics have been noted in the LBGs (e.g. Pettini et al. 2000; Pettini et al. 2001; see also Schaye 2001) and can be explained in terms of radiative transfer effects of Ly$`\alpha `$ photons escaping through outflowing gas. Photons at the Ly$`\alpha `$ wavelength in the rest-frame of the moving HI are scattered repeatedly and are eventually absorbed by dust; those in the wings of the Ly$`\alpha `$ line have a smaller scattering cross-section and escape from the neutral gas (Kunth et al. 1998). In terms of Ly$`\alpha `$ line emission, the properties of all three DLA galaxies are therefore in good agreement with those of the LBGs.

## 5 Discussion and conclusions

In this paper we have presented new STIS imaging data and VLT spectroscopic data on a small set of three DLA galaxies. Two of those were previously known, one is a new discovery. The small set of three objects does not allow a conclusive answer to the dispute over whether DLA absorbers arise in large disks or merging sub-clumps. Our data set does allow us, however, to address the more fundamental (and model independent) question of the possible relationship between DLA and Lyman-break galaxies.

We have compared the emission properties of the three DLA galaxies with the emission properties of the population of LBGs. If DLA galaxies are LBGs then the measured values for an emission property for any DLA galaxy will lie within the range of values measured for LBGs of similar absolute magnitude – regardless of the selection biases that are inherent in the sample of DLA galaxies. We have tested this prediction, making the comparison for the following emission properties: half-light radius, radial profile (Sersic $`n`$ parameter), optical-to-near-infrared colour, morphology, Ly$`\alpha `$ emission EW, and Ly$`\alpha `$ emission velocity structure. In all cases the measured values for the DLA galaxies lie within the range measured for LBGs. None of the measurements is in conflict with the prediction. Therefore we conclude that the measured emission properties of the three damped Ly$`\alpha `$ galaxies studied here are consistent with the conjecture that damped Ly$`\alpha `$ galaxies are Lyman-break galaxies.

Is the number of detections in the follow-up spectroscopy so far completed of the 41 NICMOS candidates also consistent with this conjecture? Since two of the DLA galaxies discussed here were previously known they are ignored in this statistical analysis. We have obtained deep VLT FORS spectra, average integration time 6000s, of 11 of the 22 NICMOS candidates in 6 fields, targeting 8 DLA absorbers. We have identified and confirmed one new DLA galaxy, and we would expect to have confirmed $`(8\times 11\times F_{NIC}\times F_{Ly\alpha })/22=4\times F_{NIC}\times F_{Ly\alpha }`$ new DLA galaxies. Here the factor $`F_{NIC}`$ is the fraction of DLA galaxies detectable in the NICMOS images i.e. above the magnitude and impact-parameter detection limits of our survey (approximately $`H_{160}<25`$, $`b>0.6\mathrm{}`$ – full details provided in Warren et al., 2001) and $`F_{Ly\alpha }`$ is the fraction of DLA galaxies with detectable Ly$`\alpha `$ flux (our $`5\sigma `$ detection limit ranges from 0.2 to $`2.0\times 10^{16}`$ erg s-1 cm-2 depending on redshift, column density and impact-parameter of the DLA). If DLA galaxies are LBGs then both of the factors $`F_{NIC}`$ and $`F_{Ly\alpha }`$ are likely to be substantially less than unity (§2 and §4). Therefore the detection of only one new DLA galaxy so far is consistent with the hypothesis, but the number of detections is clearly too small to rule out other interpretations.

The combination of information on both the stellar (emission) and interstellar (absorption) components of the DLA galaxies is crucial for a complete understanding of these objects. The result presented here is a step towards combining the information from absorption-line and emission studies of high-redshift galaxies. To prove that DLA galaxies are LBGs, the test we have made is a necessary but not sufficient condition. We are building a larger sample of DLA galaxies with well-defined selection properties, from spectroscopy of candidates found in deep HST images. A definitive test will be possible with this sample. Knowing the selection biases we can determine the relation between gas radius and luminosity that explains the number of DLA galaxies detected, their luminosities, and their impact parameters.

We thank Martin Haehnelt for a helpful conversation and Max Pettini, Jason Prochaska, and Hsiao Wen Chen for helpful comments on an earlier version of this manuscript. SJW is grateful to ESO for hospitality during part of the writing of this paper. Part of the reduction of the STIS data was done at STScI, Baltimore and PM gratefully acknowledges support from the cheerful STScI support staff and in particular useful discussions with Charles Proffitt on the subject of STIS psfs. SMF is grateful to ESO and the Max-Planck Institut für Astrophysik, Garching for hospitality in the early stages of this work and the Carnegie Observatories, Pasadena, in the late stages.


text:
# Jamming transitions and avalanches in the game of Dots-and-Boxes

## I The model

As a model for the game that can be both simulated with reasonable computational effort and handled analytically, we use the following rules of behaviour for the players:

1. Occupy squares if possible: in accordance with the rules that prescribe a greedy strategy, the active player searches the board for all free edges with an adjacent square with exactly three occupied edges around it, picks one of them at random, occupies it, and continues his turn. This is repeated as often as possible.
2. Create no opportunities for the opponent: if there are no squares to be taken, the player looks for free edges whose adjacent squares both have less than two occupied edges. One of those is picked at random and occupied, and the turn ends.
3. Minimize the opponent’s score: if both prescriptions (1) and (2) yield no suitable edge, the active player has to pick the third edge around some square, thus giving his opponent a chance to take it according to (1). He checks all edges to see how many points his opponent would gain from them. He picks one of those that give away the smallest number of points, occupies it, and ends the turn.

To get rid of boundary effects, the lattice on which the game is played is assumed to have periodic boundary conditions. Furthermore, the analytical treatment is strictly valid in the limit of infinite lattice size only.

The lattice has $`N`$ squares. Accordingly, the number of edges is $`2N`$. The number of turns is counted by $`T`$. The rescaled time $`t=T/(2N)`$ runs from 0 to at most 1; however, since more than one edge per turn is occupied in the end-game, the game ends at times $`t<1`$. Two other useful quantities to describe the system are the number $`P`$ of occupied edges (or $`p=P/(2N)`$, the probability that a given edge is occupied), and the number of filled squares $`S`$ (or $`s=S/N`$, respectively).

## II The early game

As mentioned in the introduction, the early game can be treated as a special case of random sequential adsortion (RSA) Evans (1993). The basic idea of RSA is to deposit particles (atoms, dimers, or, in our case, edges) on randomly chosen sites of a surface (often on a regular lattice) unless this deposition violates restrictions posed by particles that were adsorbed before. In Dots-and-Boxes, the edges form a square lattice with a peculiar short-range three-particle repulsion.

The usual procedure to treat RSA problems analytically is to solve a set of coupled differential equations that describe the probability of encountering the various possible configurations of particles (for details, see Evans (1993)). Unfortunately, this set of ODEs is generally not closed, i.e., the equations for small configurations include probabilities of larger configuration, which in turn depend on still larger configurations. At some point, one has to neglect correlations and truncate the equations. Since the number of terms increases dramatically with increasing order of truncation, we will use the simplest approximation that still captures the interaction correctly.

This approximation characterizes each free edge by two indices – the first index $`i`$ for the number of occupied edges surrounding the adjacent square above or to the left of the considered edge, the second index $`j`$ for the occupation number of the square below or to the right (for example, Fig.1 shows all configurations with index $`21`$ around a horizontal edge). One can then count the number $`F_{ij}`$ of free edges with indices $`ij`$, and determine their density $`f_{ij}=F_{ij}/(2N)`$.

All possible configurations with the same indices are now assumed to have the same probability, and correlations beyond nearest neighbours are neglected. For example, if we consider a free edge with indices $`11`$, we have no information on the free edges surrounding its adjacent squares, except that one of their indices must be 1. Accordingly, we assume that the other index (let us say, $`j`$) follows a simple conditional probability, $`\text{Prob}(j|i)=f_{ij}/_kf_{ik}`$.

In the inititial phase of the game, there is an extensive number of edges that can be occupied without giving the opponent an opportunity to score, following prescription (2) of the strategy. These are edges from the categories 00, 10, 01, and 11. Now we can see what happens to the $`f_{ij}`$ if an edge is occupied, i.e., a time step $`dt=1/(2N)`$ elapses. The chosen edge (let us say, it has indices $`kl`$) is occupied, and $`f_{kl}`$ is decreased by $`df_{kl}=1/(2N)`$. This means that for all edges that can be chosen, the differential equation for $`f_{kl}`$ includes a loss term equal to the probability that an edge with index $`kl`$ is picked:

$$\frac{df_{kl}}{dt}=\frac{f_{kl}}{f_f}+\mathrm{}\text{ for }k,l\{0,1\},$$
(1)
where $`f_f`$ is the density of edges that can be occupied before the transition,

$$f_f=f_{00}+f_{01}+f_{10}+f_{11}.$$
(2)
The free edges in the squares next to the chosen edge must be updated: the index corresponding to the considered square is increased by one. For instance, assume that the upper/left adjacent square of the chosen edge has occupation number 0 (which happens with probability $`(f_{00}+f_{01})/f_f`$). Two of the three remaining edges around that square will then have their second index increased by one, whereas one will have its first index incremented. This leads to gain and loss terms in the ODEs:

$`f_{i0}`$ $`:`$ $`2{\displaystyle \frac{f_{i0}}{_{k=0}^2f_{k0}}}{\displaystyle \frac{f_{00}+f_{01}}{f_f}};`$
$`f_{0i}`$ $`:`$ $`{\displaystyle \frac{f_{0i}}{_{k=0}^2f_{0k}}}{\displaystyle \frac{f_{00}+f_{01}}{f_f}}`$
$`f_{i1}`$ $`:`$ $`+2{\displaystyle \frac{f_{i0}}{_{k=0}^2f_{k0}}}{\displaystyle \frac{f_{00}+f_{01}}{f_f}};`$
$`f_{1i}`$ $`:`$ $`+{\displaystyle \frac{f_{0i}}{_{k=0}^2f_{0k}}}{\displaystyle \frac{f_{00}+f_{01}}{f_f}}.`$ (3)
What follows is a rather tedious summation of terms for the possible combinations of indices.

The equations derived in this fashion can be simplified considerably by assuming the symmetry $`f_{ij}=f_{ji}`$, thus keeping only the categories with $`ij`$. With the abbreviations

$`r_0`$ $`=`$ $`(f_{00}+f_{10})/(f_{00}+f_{10}+f_{20});`$ (4)
$`r_1`$ $`=`$ $`(f_{10}+f_{11})/(f_{10}+f_{11}+f_{21}),`$ (5)
the resulting system of differential equations looks as follows:

$`{\displaystyle \frac{df_{00}}{dt}}`$ $`=`$ $`(f_{00}6f_{00}r_0)/f_f;`$
$`{\displaystyle \frac{df_{10}}{dt}}`$ $`=`$ $`(f_{10}+3f_{00}r_03f_{10}r_02f_{10}r_1)/f_f;`$
$`{\displaystyle \frac{df_{11}}{dt}}`$ $`=`$ $`(f_{11}+6f_{10}r_04f_{11}r_1)/f_f;`$
$`{\displaystyle \frac{df_{20}}{dt}}`$ $`=`$ $`(3f_{20}r_0+2f_{10}r_1)/f_f;`$
$`{\displaystyle \frac{df_{21}}{dt}}`$ $`=`$ $`(2f_{21}r_1+2f_{11}r_1+3f_{20}r_0)/f_f;`$
$`{\displaystyle \frac{df_{22}}{dt}}`$ $`=`$ $`4f_{21}r_1/f_f.`$ (6)
This can be solved numerically and compared to simulations. The agreement is very good, but becomes slightly worse close to the jamming transition – the point where $`f_f`$ becomes zero and the early game ends. The accuracy could probably be improved by including probabilies of larger configurations; however, it is not worth the effort, since the calculated numbers agree with results from real play only in the order of magnitude anyway (see Section V). The predictions and numerical values for the jamming time $`t_J`$ and the order parameters at that time are given in Table 1:

## III The end-game

When the last edge from the categories 00, 01, 10, and 11 has been taken, prescription (2) is no longer an option. The game now alternates between prescriptions (1) and (3): Each player’s turn begins by filling the squares of the avalanche that his opponent has offered him by placing the third edge around some square (prescription (1)). When all possible squares have been taken, the active player now determines the avalanche that his opponent must take (prescription (3)). Since the length of the avalanche triggered by placing an edge is not a local property of that edge (and highly correlated to that of neighboring edges), a description by differential equations analogous to the early game makes little sense. However, since the state of the system at the transition largely determines the options of the players later on, we can make quantitative predictions about the end-game from the knowledge gained in Sec. II.

Fig. 3 shows the state of a game with $`15\times 15`$ squares at the jamming transition. The squares are segments of a tunnel if they have two occupied egdes. They represent tunnel branchings (or single defects, marked by empty circles) if they have only one edge, and tunnel crossings (or double defects, marked by full circles) if they have none.

Thus, three edges from the $`f_{21}`$-category form a single defect, whereas four $`f_{20}`$-edges make up a double defect. The density of defects can be calculated directly from the order parameters at the transition.

As an additional check, one can make sure that the total number of edges (both occupied and free) add up to $`2N`$, and that the number of squares (tunnels and defects) add up to $`N`$. This leads to the following equations:

$`t_J+2f_{20}+2f_{21}+f_{22}`$ $`=`$ $`1;`$ (7)
$`3f_{20}+(10/3)f_{21}+2f_{22}`$ $`=`$ $`1,`$ (8)
which are fulfilled for both the analytical and the experimental values.

This enables us to make some statements about the avalanches or chains of occupied squares that occur in the late game. An avalanche started in a tunnel fills the tunnel and ends at the defect edges on both sides of the tunnel. Thus, with $`4N(f_{21}+f_{20})`$ defect edges, there are at least $`N_A=2N(f_{21}+f_{20})`$ different potential avalanches at the time of the transition. (For now, we neglect avalanches in closed areas and other complications.) We can also calculate the number of tunnel segments (squares with occupation 2) from the order parameters: $`N_T=N(4/3)Nf_{21}Nf_{20}`$. This yields an average length of the tunnel segments of $`N_T/N_A4.5`$ for the values of $`f_{ij}`$ from the simulation, as given above.

Note that in an analogy to the “waiting time paradoxon”, the average avalanche length becomes larger than the mentioned value of 4.5 if avalanches are started at randomly chosen edges rather than randomly chosen tunnel segments, because longer tunnels include more edges and are thus chosen with higher probability:

Let us assume that the probability distribution $`P_{av}(l)`$ of tunnel lengths $`l`$ follows an exponential with a decay constant $`1/l^{}`$. Since $`l1`$, the normalization constant is $`\mathrm{exp}(1/l^{})1`$, and the average value is

$$l_{av}=\underset{l=1}{\overset{\mathrm{}}{}}l(e^{1/l^{}}1)e^{l/l^{}}=\frac{1}{1e^{1/l^{}}}.$$
(9)
With the mentioned value of $`l_{av}=4.5`$, one gets $`l^{}4.0`$ and $`P_{av}(l)0.284\mathrm{exp}(l/4.00)`$. Since each tunnel avalanche of $`l`$ squares length has $`l+1`$ edges where it can be started, the probability distribution of avalanche lengths averaged over free edges follows the form $`P_{ed}(l)(l+1)\mathrm{exp}(l/l^{})`$.

This agrees fairly well with simulations, as seen in Fig. 4. However, there is a preference for even avalanche lengths, which can partly be explained with the presence of closed areas. These are areas that contain no defects and are separated from the rest of the board by occupied edges. They include an even number $`4`$ of squares. The probability of an edge being in a closed area of size $`l`$ is shown in Fig. 4 (open circles). Even if it is taken into account, even avalanches are more likely, for reasons that are still unclear.

The density of defects at the transition allows for a prediction of the total number of turns in the end-game: When an avalanche is terminated at a defect, that defect is turned into a tunnel segment (if it was a single defect before) or into a single defect (if it was a double defect). Therefore the number of defect-terminated avalanches $`N_{DTA}`$ is half the number of single defects plus the number of double defects:

$$N_{DTA}=((2/3)f_{21}+f_{20})/N0.069N.$$
(10)
This means that the time difference from $`t_J`$ to the end of the game $`t_E`$ should be at least

$$t_Et_JN_{DTA}/(2N)=0.034.$$
(11)
However, this is really only a lower bound on the time found in simulations, $`t_Et_J0.054`$. The reason for the deviation is the existence of avalanches that do not change the number of defects, namely, avalanches in closed areas. These areas are not necessarily present at the beginning of the end-game. Instead, they may initially be half-closed areas: areas that are not quite closed, but connected to the rest of the system by a single defect. Depending on whether the avalanche is started in the tunnel outside the half-closed area or inside it, it is either turned into a closed area, or it is filled, the defect is removed, and the adjacent tunnel is filled as well. Since it is usually desirable to give the opponent as few points as possible, most of the half-closed areas in real play will be turned into closed areas, and then filled.

Apart from these exceptions, avalanches tend to get dramatically longer as the end-game goes on (see Fig. 5). This is due to two effects: first, small avalanches are triggered earlier than larger ones due to prescription (3), and thus removed; second, avalanches that stop at a single defect turn it into a tunnel segment, merging two potential avalanches into one.

## IV Other lattices

The game can be played on lattices other than the square, as long as there is a notion of an edge separating two cells, and placing a single edge is not enough to make scoring possible. The simplest case is the triangular lattice, where there are only single defects, and the relevant order parameters in the early game are $`f_{00}`$, $`f_{10}`$, and $`f_{11}`$. The corresponding differential equations are

$`{\displaystyle \frac{df_{00}}{dt}}`$ $`=`$ $`14{\displaystyle \frac{f_{00}}{f_{00}+f_{10}}};`$ (12)
$`{\displaystyle \frac{df_{10}}{dt}}`$ $`=`$ $`2{\displaystyle \frac{f_{00}f_{10}}{f_{00}+f_{10}}};`$ (13)
$`{\displaystyle \frac{df_{11}}{dt}}`$ $`=`$ $`4{\displaystyle \frac{f_{10}}{f_{00}+f_{10}}}.`$ (14)
They can even be solved analytically by introducing a rescaled time $`\tau `$ with $`d\tau =(f_{00}+f_{10})^1dt`$, and solving the remaining system of linear ODEs in $`\tau `$ with constant coefficients. One gets

$`f_{00}(\tau )`$ $`=`$ $`e^{4\tau }(2e^\tau );`$ (15)
$`f_{10}(\tau )`$ $`=`$ $`2e^{4\tau }(e^\tau 1);`$ (16)
$`f_{11}(\tau )`$ $`=`$ $`2e^{4\tau }(8/3)e^{3\tau }+2/3.`$ (17)
Using Eqs. (15) and (16), $`t`$ can be calculated:

$$t(\tau )=(1e^{3\tau })/3;\tau (t)=(1/3)\mathrm{ln}(13t).$$
(18)
Again, agreement between theory and simulation is very good, as seen in Fig. 6. The predicted and observed values for the jamming transition and the order parameters are given in Table 2.

Other possible lattices include hexagonal and three-dimensional cubic lattices. In the latter case, edges correspond to faces of unit cubes. Interestingly, the differential equations are the same for both the hexagonal and the 3D-cubic lattice, since both have six edges/faces surrounding each hexagon/cube. Although structurally simple, the equations involve fifteen order parameters and are not written out for the sake of brevity.

Of course, the whole range from single to quadruple defects can occur; however, multiple defects are rarer than single ones, as seen in Table 3.

Results from simulations confirm the picture predicted by calculations, with the usual deviations on the order of $`10^3`$.

Is the game still interesting on other lattices? Disregarding the practical difficulties of playing on a 3D-cubic lattice, all basic mechanisms of the game still work, including closed and half-closed areas. A rough estimate shows that the initial avalanche length (averaged over possible avalanches) is 4.8 for the triangular lattice and 3.6 for the 3D-cubic and hexagonal lattice, similar to the square lattice. We therefore expect that real-life games on other lattices would not be much different from Dots-and-Boxes on regular square lattices.

## V Comparison to real play

We let some coworkers play a computer version of Dots-and-Boxes (with periodic boundary conditions and $`N=10\times 10`$) to see if their style of play is well described by the assumptions in Section I. Generally speaking, human players did not place edges at random in the early game; instead they tended to add edges to existing structures. In some cases, this led to a significantly lower number of defects, and thus longer avalanches. One pair of players chose to get rid of the periodic boundary conditions by drawing a frame around the board early in the game.

Nevertheless, some games showed quantitative similarities to our theoretical predictions. The order parameters from one of these games is shown in Figs 7 and 8.

Our test players did not try tactical subtleties like adjusting the number of avalanches in order to get the last (and presumably longest) one. They were usually happy if they avoided blatant mistakes.

Of course, one could include human tendencies and extend the model to include cooperative sequential adsorption Evans (1993), where edges are preferrably placed next to edges occupied before. However, since this could not describe all human players with the same set of parameters and would probably give no qualitative new insights, the usefulness of this extension is questionable.

## VI Summary and conclusion

We gave a statistical treatment of the game of Dots-and-Boxes, using some simplifying assumptions for the behaviour of the players. In the early game, since a finite fraction of edges can be chosen, a mean-field description given by a system of coupled differential equations works well. It makes predictions about the point where avalanches start and the degree of geometrical frustration at that point. The same scheme works for all kinds of regular lattices; the relevant quantity is the number of edges or faces around a cell, such that hexagonal and 3D-cubic lattices are described the same equations.

These predictions allow for statements regarding the statistics of avalanches, as well as the total number of turns in the end-game. The presence of closed and half-closed areas makes the situation more complicated; unfortunately, they cannot be captured by the approximations made in the calculation of the early game.

While results from calculations give good agreement with simulations, human players have various habits that cannot be easily included in an all-encompassing mean-field treatment (“I like making corners. They look nice.”) Thus, while our analysis has yielded some insight in the underlying processes of Dots-and-Boxes, quantitative agreement with human play is not always satisfactory.

###### Acknowledgements.

R. M. acknowledges finacial support by the GIF. We also like to thank our test players, and Wolfgang Kinzel and Michael Biehl for pointing us to relevant literature.