text:
# Entanglement between a qubit and the environment in the spin-boson model

## Abstract

The quantitative description of the quantum entanglement between a qubit and its environment is considered. Specifically, for the ground state of the spin-boson model, the entropy of entanglement of the spin is calculated as a function of $`\alpha `$, the strength of the ohmic coupling to the environment, and $`\epsilon `$, the level asymmetry. This is done by a numerical renormalization group treatment of the related anisotropic Kondo model. For $`\epsilon =0`$, the entanglement increases monotonically with $`\alpha `$, until it becomes maximal for $`\alpha 1^{}`$. For fixed $`\epsilon >0`$, the entanglement is a maximum as a function of $`\alpha `$ for a value, $`\alpha =\alpha _M<1`$.

Due to the promise of quantum computation there is currently considerable interest in the relationship between entanglement, decoherence, entropy, and measurement. Motivated by quantum information theory several authors have recently investigated entanglement in quantum many-body systems osborne ; nature ; hines ; vidal . It is often stated that decoherence or a measurement causes a system to become entangled with its environment. The purpose of this paper is to make these ideas quantitative by a study of the simplest possible model, the spin-boson model leggett.87 ; weiss . This describes a qubit (two-level system) interacting with an infinite collection of harmonic oscillators that model the environment responsible for decoherence and dissipation. Specifically, we show how the entanglement between a superposition state of the qubit and the environment changes as the coupling between the qubit and environment increases. One interesting result is that we find that the qubit becomes maximally entangled with the environment when the coupling $`\alpha `$ approaches a particular finite value ($`\alpha 1^{}`$). Furthermore, at this value the model undergoes a quantum phase transition, which is consistent with recent observations that often entanglement is largest near quantum critical pointsosborne ; nature ; hines ; vidal .

The spin-boson model. The Hamiltonian is leggett.87 ; weiss

$`H_{SB}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }\sigma _x+{\displaystyle \frac{1}{2}}\epsilon \sigma _z+{\displaystyle \underset{i}{}}\omega _i(a_i^{}a_i+{\displaystyle \frac{1}{2}})`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}\sigma _z{\displaystyle \underset{i}{}}\lambda _i(a_i+a_i^{}),`$
where $`\mathrm{\Delta }`$ is the bare tunneling amplitude between the two quantum mechanical states $``$ and $``$, $`\epsilon `$ is the level asymmetry (or bias), $`\omega _i`$ are the frequencies of the oscillators and $`\lambda _i`$ the strength with which they couple to the two quantum mechanical states. The effect of the oscillator bath is completely determined by the spectral function $`J(\omega )`$, defined belowleggett.87 . We will only consider the ohmic case, where it is has a linear dependence on frequency

$$J(\omega )=\pi \underset{i}{}\lambda _i^2\delta (\omega \omega _i)=2\pi \alpha \omega $$

for $`\omega \omega _c`$, and $`\alpha `$ is the dimensionless dissipation strength. The cutoff frequency, $`\omega _c\mathrm{\Delta }`$. This model can describe the decoherence of Josephson junction qubits, such as those recently realized experimentallyexpt , due to voltage fluctuations in the electronic circuitschon , and $`\alpha `$ can be expressed in terms of resistances and capacitances in the circuit and so this is an experimentally tunable parameter. Recent results show it is possible to construct devices, with $`\alpha 1`$, the regime required for quantum computation. However, when modeling measurements one has $`\alpha 1`$.

The dynamical properties of the model have been extensively studied. In particular, suppose the spin (qubit) is initially in a pure state which is a product state of up spin and the environment state, then the coherent Rabi oscillations that would be observed in the absence of coupling to the environment are modified as follows. One finds distinct behaviour for $`0<\alpha <1/2`$ (damped coherent oscillations), $`1/2<\alpha <1`$ (exponential decay), and $`1<\alpha `$ (localization; i.e., the spin remains in the up state) leggett.87 ; weiss ; saleur ; costi.96 .

Entropy of entanglement. We now consider a quantitative description of the entanglement of the qubit with the environment. A good entanglement measure for a pure state is the von Neumann entropy or entropy of entanglementbennett ; vedral

$$E(\rho )=Tr(\rho \mathrm{log}_2\rho )$$
(2)
where $`\rho `$ is the reduced density matrix of the qubit. This is a two by two matrix given by

$$\rho =\frac{1}{2}\left(1+\underset{a=x,y,z}{}\sigma _a\sigma _a\right),$$
(3)
where $`\sigma _a`$ denotes the expectation value in the state of interest. In this case Eq. (2) reduces to

$$E(\rho )=p_+\mathrm{log}_2p_+p_{}\mathrm{log}_2p_{},$$
(4)
where $`p_\pm `$ are the eigenvalues of the density matrix,

$$p_\pm =\frac{1}{2}\left(1\pm |\stackrel{}{\sigma }|\right).$$
(5)
For $`\epsilon =0`$ the only nonzero value of $`\sigma _a`$ is $`\sigma _x`$. At $`T=0`$ it is given by

$$\sigma _x=2\frac{E_0}{\mathrm{\Delta }},$$
(6)
where $`E_0`$ is the ground-state energy of the Hamiltonian (1), and use has been made of the Feynman-Hellmann theorem. That the other values are zero can be seen by symmetry as follows. In general the Hamiltonian is invariant under the reflection in spin space, $`\sigma _y\sigma _y`$. Hence, all eigenstates must have a definite parity under this transformation. Thus, $`\sigma _y=\sigma _y`$ for all states and so $`\sigma _y=0`$ at any temperature. For $`\epsilon =0`$ the Hamiltonian is also invariant under the joint transformation $`\sigma _z\sigma _z`$ and $`a_ia_i`$ and so $`\sigma _z=0`$ at any temperature, provided there is no symmetry breaking.

The challenge is now to evaluate the ground-state expectation values $`\sigma _x`$ and $`\sigma _z`$. For $`\alpha >1/2`$ and particularly for $`\alpha 1`$ this is a highly nontrivial problem because in this regime nonperturbative effects become importantleggett.87 ; weiss . However, we show how these expectation values can be evaluated using the numerical renormalization group (NRG) applied to the equivalent anisotropic Kondo model.

Anisotropic Kondo Model. The above model is equivalent to the anisotropic Kondo model (AKM), defined by anderson.69

$`H`$ $`=`$ $`{\displaystyle \underset{k,\sigma }{}}ϵ_kc_{k\sigma }^{}c_{k\sigma }+{\displaystyle \frac{J_{}}{2}}(c_0^{}c_0S^{}+c_0^{}c_0S^+)`$ (7)
$`+`$ $`{\displaystyle \frac{J_{}}{2}}(c_0^{}c_0c_0^{}c_0)S^z+g\mu _BhS_z.`$
The first term represents a free electron conduction band. We use a flat density of states $`\rho _0=1/2D_0`$ per spin, with $`2D_0`$ the bandwidth. The second and third terms represent the transverse and longitudinal parts of the exchange interaction between a $`S=1/2`$ impurity and the local conduction-electron spin density, and the last term represents a Zeeman term for a magnetic field coupling only to the impurity spin. The correspondence between $`H`$ and $`H_{SB}`$, established via bosonizationguinea.85 ; costi.99 , implies $`\epsilon =g\mu _Bh`$, $`\frac{\mathrm{\Delta }}{\omega _c}=\rho _0J_{}`$ and $`\alpha =(1+2\delta /\pi )^2`$, where $`\mathrm{tan}\delta =\pi \rho _0J_{}/4`$. $`\delta `$ is the phase shift for scattering of electrons from a potential $`J_{}/4`$ leggett.87 ; costi.96 ; guinea.85 . We choose $`\omega _c=2D_0`$ so that $`\mathrm{\Delta }=J_{}`$. This equivalence has been used extensively to make predictions about the dynamicscosti.96  and thermodynamics costi.99  of the ohmic spin-boson model. The relevant low energy scale for the thermodynamics is the Kondo scale $`T_K(J_{},J_{})`$ which is identified with the renormalized tunneling amplitude, $`\mathrm{\Delta }_r`$, of the spin-boson model,

$$\frac{\mathrm{\Delta }_r}{\omega _c}=\left(\frac{\mathrm{\Delta }}{\omega _c}\right)^{1/(1\alpha )}.$$
(8)
We restrict ourselves in this paper to the longitudinal sector of the AKM, i.e., to $`J_{}<|J_{}|`$, where the simple parameter correspondence between the models given above remains valid to lowest order in $`\mathrm{\Delta }/\omega _c=\rho _0J_{}`$. For larger values of $`\mathrm{\Delta }/\omega _c`$, $`\alpha `$ will acquire a renormalization due to finite $`\mathrm{\Delta }=J_{}`$, as indicated by the scaling analysis of the AKM in Refs. anderson.69 ; costi.99 . This renormalization, however, is important mainly for the transverse sector of the AKM, $`J_{}>|J_{}|`$, which we do not consider in this paper.

We turn now to the evaluation of $`\sigma _x`$. The equivalence between models ensures that the AKM has (to within an additive constant) the same ground-state energy $`E_0`$ as that of the spin-boson model. At $`T=0`$, we therefore find, in analogy to Eq. (6) applied to the AKM with $`\mathrm{\Delta }=J_{}`$, that

$$\sigma _x=c_{0,}^{}c_{0,}S^{}+H.c.$$
(9)
i.e., $`\sigma _x`$ can be obtained from a local static correlation function. Another way of seeing that this relation is valid, is to note that the unitary transformation in bosonization which transforms $`H`$ into $`H_{SB}`$ also transforms $`(c_{0,}^{}c_{0,}S^{}+H.c.)`$ into $`\sigma _x`$ of the spin-boson model (details of this mapping can be found in Ref. guinea.85  and in greater detail in Appendix A of Ref. costi.99, ). The same unitary transformation on the AKM transforms $`S_z`$ into $`\sigma _z/2`$ of the spin-boson model. The latter can therefore be calculated directly within the AKM as a thermodynamic average $`2S_z`$.

Method. The above local correlation function can be calculated from Wilson’s NRG method wilson.75+kww.80  which has been shown to give very reliable results for quantum impurity models such as the AKMdmrg . The approach used here allows in addition the calculation of local dynamical quantities, such as the dynamical susceptibility $`c_{0,}^{}c_{0,}S^{};c_{0,}^{}c_{0,}S^+`$ costi.94 . In outline (see Ref. wilson.75+kww.80  for the details), the procedure consists of introducing a logarithmic mesh of $`k`$ points $`k_n=\mathrm{\Lambda }^n,\mathrm{\Lambda }>1`$ for the conduction band and performing a unitary transformation of the $`c_{k\sigma }`$ such that $`f_{0\sigma }=_kc_{k\sigma }`$ is the first operator in a new basis, $`f_{n\sigma },n=0,1,\mathrm{}`$, which tridiagonalizes $`H_c=_{k\mu }ϵ_{k\mu }c_{k\mu }^{}c_{k\mu }`$ in k space. The Hamiltonian (7) with the new discretized form of the kinetic energy is now diagonalized by the following iterative process: (a) one defines a sequence of finite size Hamiltonians $`H_N=_\mu _{n=0}^{N1}\xi _n\mathrm{\Lambda }^{n/2}(f_{n+1\mu }^{}f_{n\mu }+H.c.)+\frac{J_{}}{2}(f_0^{}f_0S^{}+f_0^{}f_0S^+)+\frac{J_{}}{2}(f_0^{}f_0f_0^{}f_0)S^z`$ for $`N0`$ and $`\xi _n1`$ for $`n1`$ wilson.75+kww.80 ; (b) the sequence of Hamiltonians $`H_N`$ for $`N=0,1,\mathrm{}`$ is iteratively diagonalized within a product basis of, typically, up to 1200 states for each iteration, up to a maximum value $`N=N_m`$. This gives the excitations and many body eigenstates at a corresponding set of energy scales $`\omega _N`$ defined by the lowest scale $`\omega _N=\mathrm{\Lambda }^{\frac{N1}{2}}`$ in $`H_N`$. The matrix elements $`m|O_{x,z}|n_N`$ for the operators $`O_x=c_{0,}^{}c_{0,}S^{}`$ and $`O_z=S_z`$, required to calculate $`\sigma _x`$ and $`\sigma _z`$, are also calculated iteratively. The choice of $`N_m`$ depends on the Kondo scale $`T_K=\mathrm{\Delta }_r`$ and hence on $`\alpha `$, but for given $`\alpha `$ (i.e. for given $`J_{},J_{}`$) should be large enough such that $`\omega _{N_m}\mathrm{\Delta }_r`$. A discretization parameter $`\mathrm{\Lambda }=1.5`$ was used throughout and we checked that the above expectation values remained unchanged on further increasing $`N_m`$. This suggested that our approximation of using a finite $`N_m`$ to calculate the thermodynamic expectation values of the infinite system is a very good one.

Symmetric case. For $`\epsilon =0`$ and $`\alpha <1`$ only $`\sigma _x`$ is nonvanishing. We show this in Fig. (1) at $`T=0`$ versus dimensionless dissipation strength $`\alpha `$ in the range $`0<\alpha <1`$, and for several values of the dimensionless tunneling amplitude $`\mathrm{\Delta }/\omega _c`$. The limiting noninteracting value $`\sigma _x1`$ is recovered as $`\alpha 0`$. In the limit $`\mathrm{\Delta }/\omega _c0`$ it vanishes at the quantum critical point of the spin-boson model $`\alpha =1`$ where $`\mathrm{\Delta }_r0`$. For any finite fixed $`\mathrm{\Delta }/\omega _c`$, however, our use of the AKM implies that the critical behaviour occurs at $`\alpha _c>1`$ with $`\alpha _c1`$ as $`J_{}0`$ (specifically, this critical behaviour occurs at the ferromagnetic-antiferromagnetic boundary $`J_{}=J_{}`$). Figure (1) also shows the $`T=0`$ entropy of entanglement of the qubit. The entropy vanishes as $`\alpha 0`$ and approaches its maximum value as $`\alpha 1^{}`$ (see also Ref. sidles.03  for weak dissipation results for E). For $`\alpha >1`$, we are in the ferromagnetic sector of the AKM where $`\sigma _z=1,\sigma _x=0`$, and the reduced density matrix eigenvalues $`p_\pm =0,1`$ giving $`E=0`$, i.e., $`E(\alpha )`$ drops discontinuously at the quantum critical point $`\alpha =1`$ long-paper . It is interesting that for a spin qubit coupled to two bosonic baths it is possible to remain in the delocalized phase (i.e., $`\sigma _z=0`$) for all dissipation strengths castro-neto.03 . Finally, we note that the entropy of entanglement is quite different from the thermodynamic entropy of the boundary (or impurity spin entropy). The latter is usually defined as $`S(\alpha )S(\alpha =0)`$ where $`S(\alpha )`$ is the total thermodynamic entropy of the system costi.99 . The impurity spin entropy is zero for $`\alpha <1`$ because the ground state of the AKM is a spin singlet for $`J_{}>0`$.

Asymmetric case. For $`\epsilon >0`$, $`\sigma _z`$ acquires a finite value analogous to the magnetization $`S_z`$ in a local magnetic field $`g\mu _Bh=\epsilon `$ in the AKM. The entanglement entropy $`E`$ now depends on $`|\stackrel{}{\sigma }|=(|\sigma _x|^2+|\sigma _z|)^{1/2}`$ via Eq. (3) and is shown in Fig. (2). The behaviour of $`E`$ as a function of $`\alpha `$ and $`\epsilon `$ is understood from the behaviour of $`\sigma _x`$ and $`\sigma _z`$ shown in Fig. (3). In particular, we now find that for arbitrary small $`\epsilon `$, the entanglement entropy first increases with increasing $`\alpha `$ before reaching a maximum value at $`\alpha =\alpha _M<1`$ and then decreasing as $`\alpha 1`$. This behaviour arises from the competition between the $`\alpha `$ dependence of $`\sigma _x`$ and $`\sigma _z`$ in Fig. (3). Whereas $`\sigma _x`$ continues to decrease monotonically with increasing $`\alpha `$ (as for $`\epsilon =0`$), it is seen that $`\sigma _z`$ increases monotonically with increasing $`\alpha `$ with $`\sigma _z1`$ as $`\alpha 1`$. The condition for full polarization, $`\sigma _z1`$, is $`\epsilon \mathrm{\Delta }_r`$. For any $`\epsilon >0`$, this condition is always satisfied since $`\mathrm{\Delta }_r0`$ as $`\alpha 1`$. It follows that $`|\stackrel{}{\sigma }|`$ has a minimum as a function of $`\alpha `$ and that the entanglement entropy, for $`\epsilon >0`$, has a maximum at $`\alpha =\alpha _M<1`$ before decreasing again as $`\alpha 1`$.

Finally, we suggest several directions for future work.

(i) This work focused solely on static properties of the spin-boson model. It would be interesting to consider dynamics, for example, the longitudinal and the transverse dynamical susceptibilities, and hence extract the decoherence and relaxation rates for an ohmically coupled qubit. In addition, it is interesting to ask how the entanglement varies with time if the initial state has no entanglement of the qubit and the environment.

(ii) The AKM is integrable by the Bethe ansatztsvelick.83 . The AKM can also be related to a free boson field theory with a boundary sine-Gordon termsaleur ; saleur2  which is also integrable by the Bethe ansatz. Exact expressions can be obtained for the free energy. It involves solving a set of thermodynamic Bethe ansatz (TBA) equations. At $`T=0`$ the impurity ground state energy is going to be related to $`T_K`$. The real problem is getting results for arbitrary anisotropies (dissipation strengths)TBA .

(iii) Recently it was shownvidal  that if in a quantum critical system one calculates the entropy of entanglement of a subsystem of size $`L`$ with the rest of the system this equals the geometric entropy previously calculated for the corresponding conformal field theory (motivated by questions concerning black hole thermodynamics!)holzhey . It would be interesting to perform similar calculations for the relevant boundary field theory.

(iv) The NRG can also be used to reliably calculate properties of the spin-boson model at nonzero temperaturelong-paper . However, calculating the entanglement at nonzero temperature is an open problem because it involves a mixed state and it is not practical to evaluate the measure of entanglement that has been proposed for such statesbennett .

T.A.C. acknowledges support from the Deutsche Forschungsgemeinschaft through the Sonderforschungsbereich 195. Work at UQ was supported by the Australian Research Council. We thank A. Hines, A. W. W. Ludwig, G. J. Milburn, and H. -Q. Zhou for very helpful discussions.


format: text

text:
# Untitled Document

COMMISSIONS 27 AND 42 OF THE IAU

INFORMATION BULLETIN ON VARIABLE STARS

Number 5xxx

Konkoly Observatory Budapest 00 Month 200x HU ISSN 0374 – 0676

SPECTROSCOPIC BINARIES IN THE OPEN CLUSTER TRUMPLER 16 REVISITED

LUNA, G.J.1; LEVATO, H.2; MALARODA, S.2; GROSSO, M.2

1 Instituto de Astronomia e Geofisica, USP, Rua de Matão 1226, São Paulo, Brazil, e-mail: gjmluna@astro.iag.usp.br 2 Complejo Astronómico El Leoncito, Av.España 1512 Sur, San Juan, Argentina

The open cluster Trumpler 16 is near the Carina nebulae (NGC 3372). There is a large amount of literature on the clusters in this region. The spectral morphology of members of Trumpler 16 has been studied by Levato & Malaroda (1982) and references therein. The cluster distance has been determined using different methods, and there is an agreement around 3.5 Kpc. This value is obviously dependeds on the interstellar absorption which may be abnormal for this cluster. Trumpler 16 is a rich cluster with a large variety of objects like WR stars, O3 stars, spectroscopic binaries and the singular $`\eta `$ Carinae. The age estimated for Trumpler 16 is aproximately 107 years (Massey & Johnson, 1993).

The study of spectroscopic binary systems is needed in order to tackle one of the most classic problems in astronomy namely the stellar mass determination. A good knowledge of mass for individual stars allows its calibration with other stellar parameters which are most easily observable, and also permits to test current stellar evolution theories. Spectroscopic binaries in open clusters are important in yet another context. It has been proposed that the average projected rotational velocity of the stars in an open cluster depends on its content of close spectroscopic binaries and magnetic chemically peculiar stars (Abt & Sanders, 1973). Both tidal and magnetic braking are responsible for reducing the axial rotation. (Levato et al, 1987; Abt et al. 1973)

We used the Jorge Sahade-2.1m Telescope at CASLEO, San Juan, Argentina. A REOSC echelle spectrograph was employed during the following nights: February 27th 1996, from April 21st to 23rd 2000 and from March 10th to 11th 2001. The spectra centered on the blue wavelength region were recorded on a TEK1024 CCD and the one pixel resolution is 0.14 Å. The usual flat-fields and bias frames were obtained each night and the wavelength calibration was done using a Th-Ar lamp. The reduction was made with the standard procedure, using IRAF111IRAF is distributed by the National Optical Astronomy Observatories, operated by AURA, Inc., under cooperative agreement with the NSF.

We have selected five stars for this project, probable members of Trumpler 16, from Levato et al. (1991). Three of them were chosen because the number of observations in their radial velocity curves was smaller than ten, and the other two because the authors considered them as radial velocity variables. To derive the radial velocities we have measured the Doppler shifts of H, He I, He II, Si IV, C IV, N III, C III lines. (Walborn et al. 2000), and applied an heliocentric correction to the measurements. The observational results for the three stars for which we have recalculated the orbits can be seen the Table 1, where we have indicated in successive columns the Julian Date, the average radial velocity for each spectra, the number of lines measured, and the probable errors. The stars are identified by their HD number or Feinstein, Marraco, and Muzzio (1973) numbers.

Table 1: Radial Velocity Observations

| Identification | Julian Date | RV | n | P.E. |
| --- | --- | --- | --- | --- |
|  |  | (Km/s) | n | (Km/s) |
|  | 2450141.73 | -85.39 | 6 | 2.0 |
|  | 2451656.52 | 18.46 | 7 | 6.2 |
| #112 | 2451657.52 | 54.45 | 8 | 5.8 |
|  | 2451657.61 | 50.24 | 7 | 7.0 |
|  | 2451658.57 | -30.15 | 8 | 10.7 |
|  | 2451979.76 | -76.7 | 9 | 11.5 |
|  | 2451979.78 | -118.91 | 4 | 15.5 |
| HD93161 | 2451979.79 | -91.11 | 12 | 8.6 |
|  | 2451980.61 | -105.69 | 6 | 8.8 |
|  | 2451980.65 | -96.64 | 8 | 12.3 |
|  | 2451980.77 | -97.58 | 10 | 9.8 |
|  | 2451979.74 | -181.2 | 7 | 6.1 |
|  | 2451979.77 | -174.7 | 7 | 12.2 |
| #104 | 2451980.57 | -19.7 | 7 | 5.8 |
|  | 2451980.63 | 2.3 | 5 | 10.3 |
|  | 2451980.74 | 109.2 | 7 | 18.6 |
|  | 2451980.79 | 123.4 | 6 | 25.6 |

We have added the new observations presented in Table 1 to those published by Levato et al (1991) and searched for periods, using Morbey’s code (Morbey, 1978) and, when successful, computed orbital elements starting from the results from Levato et al. (1991) improving them with the code of Bertiau & Grobben (1969).

For HD93161, #112 and #114 stars, we found new orbital parameters, while for #10 and #110, we could not find any significant evidence of variability. We have applied an analysis of variance test (Conti et al. 1977) which confirms that the distribution of the observations for stars #110 and #10 do not depart significantly from a random one. The new orbital elements for stars HD 93161, #112, and #114 are shown in Table 2.

Table 2: Orbital Elements obtained in this work.

| Element | HD93161 | #104 | #112 |
| --- | --- | --- | --- |
| a sin i (Km) | 3.72 x 106 | 4.03 x 106 | 3.67 x 106 |
| K (Km/s) | 50.9$`\pm `$5.1 | 162.8 $`\pm `$ 11.4 | 86.0 $`\pm `$ 4 |
| e | 0.309 $`\pm `$0.116 | 0.156 $`\pm `$0.057 | 0.249 $`\pm `$ 0.026 |
| $`\omega ({}_{}{}^{})`$ | 184.4 $`\pm `$ 17.8 | 85.7 $`\pm `$ 28.4 | 286.6 $`\pm `$8.8 |
| T0(J.D.)(2.400.000+) | 45778.39 $`\pm `$ 0.22 | 45779.37 $`\pm `$ 0.13 | 45773.02 $`\pm `$ 0.1 |
| P(days) | 5.60486 $`\pm `$ 9x10-5 | 1.82303 $`\pm `$1x10-5 | 4.07997 $`\pm `$2x10-5 |
| V0 (Km/s) | -44.5 $`\pm `$ 3.4 | -32.7 $`\pm `$ 7.5 | -19.1 $`\pm `$ 1.7 |
| Mass Function | 0.06 $`\pm `$ 0.03 | 0.78 $`\pm `$ 0.18 | 0.24 $`\pm `$ 0.02 |

Figures 1, 2 and 3 show the radial velocity curves for #104, #112 and HD93161 respectively.

We have recalculated new orbital elements for three spectroscopic binary systems which belong to the open cluster Trumpler 16 and obtained the mass function for these three systems. Also we could not find significant radial velocity variations for #110 and #10 of the same cluster.

Support from Universidad Nacional de San Juan and CASLEO is deeply appreciated.

References:

Abt H.& Sanders W., 1973, ApJ, 186, 177

Bertiau F. C. & Grobben J., 1969, Specola Astr. Vaticana Ric. Astr., 8, No.1

Conti P., Garmani C. & Hutchings J., 1977, ApJ, 215, 561

Feinstein A., Marraco H. G. & Muzzio J. C., 1973, A&A, 12, 331

Levato H., Malaroda S., 1982, PASP, 94, 807

Levato H., Malaroda S., Morrell N., Garcia B., Hernández C., 1991, ApJS, 75, 869

Levato H., Malaroda S., Morrell N. & Solivella G., 1987, ApJS, 64, 487

Massey, P. & Johnson, J., 1993, Astron. J., 105, 980

Morbey, Ch., 1978, Pub. Dom. Ap. Obs., 15, 105

Walborn N. R., Lennon D., Heap S., Lindler D. Smith L., Evans J., Parker J., 2000, STSI Preprint Series, 1437


text:
# Hyperspherical partial wave calculation for double photoionization of the helium atom at 20 eV excess energy

## 1 Introduction

There has been a significant development in the last one decade in the theoretical study of double photoionization (DPI) of the helium atom. This was possible because of rapid developments in the experimental side by several groups, extending over several countries. Double photoionization of the helium atom is one of the most basic atomic processes. Even then, this problem merits further detailed studies as it involves complex three body effects and electron correlations that are not yet fully understood. Total DPI cross section for helium has been measured by several groups  and also calculated theoretically by many others . In general there is good agreement between the observed and the calculated cross sections. However, there are also some discrepancies. For example experimental results of Samson et aland of Bizau and Wuilleumier, although agree within 15 percent of each other, there are systematic differences at energies above 1 Rydberg. So also is the case regarding the theories. Results of hyperspherical $``$-matrix with semi-classical outgoing waves (H$``$M-SOW) of Selles et al agree favourably with the results of Bizau and Wuilleumier  while the 2SC calculation of Pont and Shakeshaft  favours the measurements of Samson *et al*  (see ref. ). As regards the differential cross-sections, understanding of the results has much improved with the availability of results of some elaborate calculations like the time-dependent close coupling (TDCC) calculation of Colgan *et al* , the H$``$M-SOW calculation of Malegat et aland Selles *et al*  in addition to the earlier results () and the recent results ( of CCC theory, the most extensively applied theory to the problem to date. For the status of the CCC and some other methods the review article by Briggs and Schimdt  may be seen. Even then there remains the task of improving the theoretical estimates, particularly at low energies and for unequal-energy-sharing geometries (for which the results appear comparatively less satisfactory) as well as for better understanding of the physics of the DPI problem of helium

In the time-independent frame work the solution of the problem depends basically on the accurate computation of the T-matrix element given by

$$T_{fi}=\mathrm{\Psi }_f^{()}|V|\mathrm{\Phi }_i,$$
(1)
where $`\mathrm{\Phi }_i(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ is the helium ground state wave function, V is the interaction term given by

$$V=\stackrel{}{ϵ}\stackrel{}{D}.$$
(2)
$`\stackrel{}{D}`$ is the dipole operator and is given $`\stackrel{}{D}=\stackrel{}{}_1+\stackrel{}{}_2`$ (velocity form) or $`\omega _i(\stackrel{}{r}_1+\stackrel{}{r}_2)`$ (length form) and $`\mathrm{\Psi }_f^{()}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ is the final channel continuum wave function with incoming wave boundary condition for the two outgoing electrons and $`\omega _i`$ is the incident photon energy. Here $`\stackrel{}{ϵ}`$ is the photon polarization direction and $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$ are the co-ordinates of the two outgoing electrons, the nucleus being at the origin.

For accurate cross section results one needs accurate wave functions $`\mathrm{\Phi }_i`$ and $`\mathrm{\Psi }_f^{()}`$. Accurate bound state helium wave functions are easily available. There exists a number of such wave functions for the ground state (and low-lying excited states) in analytic form of different accuracies, such as simple Hartree Fock type wave function used by Maulbetsch and Briggs  or a Hylleraas type wave function given by Chandrashekhar and Herzberg  or by Hart and Herzberg . If needed, one may also readily generate an arbitrarily accurate bound state wave function along the line developed by Pekeris . But for $`\mathrm{\Psi }_f^{()}`$ there are no such simple accurate wavefunctions in analytical form. Most wave functions used in the literature are either incorrect at finite distances or in the asymptotic domain, as 3C  or 2SC  wave functions. There are now many calculations of varied accuracies depending mainly on the use of different final channel wave functions.

Without caring for the explicit form of the wave functions, Huetz and co-workers  established dependence of cross sections on various angular variables of the outgoing electrons and on energy. For the equal energy sharing geometry case it turns out to be rather simple in form. On least squares fitting this gives good representation of the triple differential cross section(TDCS) results (some authors prefer the name five fold differential cross sections (FDCS)which is more appropriate and relevant in view of some recent experimental results . However, we will continue to call it TDCS).

For the study of TDCS close to threshold there are the Wannier calculations by Faegin . These give good representation of the shape of TDCS results at 6 eV excess energy but miserably fail at higher energies.

There are also a number of detailed calculations by Maulbetsch and Briggs  which used for the final state wave function, the 3C wave function of BBK theory  and produced moderately accurate cross section results. It is well known that the 3C wave function is correct in the asymptotic domain (when all of $`r_1`$, $`r_2`$ and $`r_{12}`$ are large) but it is not accurate enough at finite distances. Similar calculations are reported by Pont and Shakeshaft . They used screened coulomb (2SC) wave functions (for the outgoing electrons) which is supposed to be a better wave function (but not asymptotically exact). The results are much better.

Later Lucey et al  tried various initial state and final state wave functions, including the 3C wave function (none accurate enough), tested gauge dependence and found much discrepancies in the results.

In this work, we consider yet another high level computational scheme which is capable of yielding reliable cross sections. This is the Hyperspherical partial wave approach of one of the authors (Das , Das et al ) which is very successful in representing the three-particle continuum wave function in the final channel very accurately.

Here we use a 20-term correlated Hylleraas type wave function of the form

$$\mathrm{\Phi }_i(r_1,r_2)=\mathrm{}e^{\frac{k}{2}(r_1+r_2)}\underset{n_1,n_2,n_3}{}C_{n_1n_2n_3}(r_1+r_2)^{n_1}(r_1r_2)^{n_2}r_{12}^{n_3}$$
(3)
given by Hart and Herzberg  for the helium ground state and use the final channel wave function calculated in the hyperspherical partial wave approach as indicated below.

It may be noted here that our present approach and the H$``$M-SOW approach are similar in some respects, although ours is a fully quantal approach whereas the H$``$M-SOW approach is partly semi classical. From the recent calculations of Selles *et al* on H$``$M-SOW approach, it is clear that for unequal energy sharing kinematics one needs to consider the asymptotic range parameter $`R_{\mathrm{}}`$ (considered below) to have values of several thousands a.u. for converged results, whereas for equal energy sharing cases only few hundreds suffice. So we restrict our present study to equal energy sharing cases only, which is much easier.

## 2 Hyperspherical Partial Wave Approach

In this section we outline the most salient features of this method. For the final state $`\mathrm{\Psi }_f^{()}`$, which needs more accurate treatment, we use hyperspherical co-coordinates $`R=\sqrt{r_{1}^{}{}_{}{}^{2}+r_{2}^{}{}_{}{}^{2}}`$, $`\alpha =atan(r_2/r_1),\widehat{r_1}=(\theta _1,\varphi _1),\widehat{r_2}=(\theta _2,\varphi _2)`$ and $`\omega =(\alpha ,\widehat{r_1},\widehat{r_2})`$ and put $`P=\sqrt{p_{1}^{}{}_{}{}^{2}+p_{2}^{}{}_{}{}^{2}}`$, $`\alpha _0=atan(p_2/p_1),\widehat{p_1}=(\theta _{p_1},\varphi _{p_1}),\widehat{p_2}=(\theta _{p_2},\varphi _{p_2})`$ and $`\omega _0=(\alpha _0,\widehat{p_1},\widehat{p_2})`$, $`\stackrel{}{p_1},\stackrel{}{p_2}`$ being momenta of the two outgoing electrons of energies $`E_1`$ and $`E_2`$ and coordinates $`\stackrel{}{r_1},\stackrel{}{r_2}`$ . We expand $`\mathrm{\Psi }_f^{()}`$ in hyperspherical harmonics (Das , Lin ) which are functions of the above five angular variables and depend on the variables $`\mathrm{}_1,\mathrm{}_2,n,L,M`$ (collectively called $`\lambda `$) which are respectively the angular momenta of the two electrons, the order of the Jacobi polynomial and the total angular momentum and its projection, in addition to the dependence on S, the total spin. It may be noted that $`L,S,\pi `$ (the parity) are conserved here.

Thus we decompose the symmetrized wave function $`\mathrm{\Psi }_{fs}^{()}`$ as

$$\mathrm{\Psi }_{fs}^{()}(R,\omega )=\sqrt{\frac{2}{\pi }}\underset{\lambda }{}\frac{F_\lambda ^s(\rho )}{\rho ^{\frac{5}{2}}}\varphi _\lambda ^s(\omega ),$$
(4)
following the expansion of the symmetrized plane wave

$`[exp(i\stackrel{}{p_1}\stackrel{}{r_1}+i\stackrel{}{p_2}\stackrel{}{r_2})`$ $`+`$ $`(1)^sexp(i\stackrel{}{p_2}\stackrel{}{r_1}+i\stackrel{}{p_1}\stackrel{}{r_2})]/(2\pi )^3`$ (5)
$`=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \underset{\lambda }{}}i^\lambda {\displaystyle \frac{j_\lambda ^s(\rho )}{\rho ^{\frac{3}{2}}}}\varphi _\lambda ^s(\omega _0)\varphi _\lambda ^s(\omega ).`$
Here $`\lambda =\mathrm{}_1+\mathrm{}_2+2n`$ and $`\rho =PR`$.

The $`F_\lambda ^s`$ satisfy an infinite coupled set of equations

$$\left[\frac{d^2}{d\rho ^2}+1\frac{\nu _\lambda (\nu _\lambda +1)}{\rho ^2}\right]F_\lambda ^s(\rho )+\underset{\lambda ^{}}{}\frac{2\alpha _{\lambda \lambda ^{}}^s}{P\rho }F_\lambda ^{}^s(\rho )=0$$
(6)
where

$$\alpha _{\lambda \lambda ^{}}^s=\varphi _\lambda ^s|C|\varphi _\lambda ^{}^s,$$
(7)
$$C=\frac{1}{cos\alpha }\frac{1}{sin\alpha }+\frac{1}{|\widehat{r_1}cos\alpha \widehat{r_2}sin\alpha |}$$
(8)
and $`\nu _\lambda =\lambda +\frac{3}{2}`$ (note that we use $`\lambda `$ with two different meanings depending on the context). Here the radial waves with L=1, S=0 and $`\pi =1`$ are relevant for the T-matrix calculations.

So here we fix $`\mu =(L,S,\pi )`$ with L=1, S=0 and $`\pi =1`$, call it $`\mu _0`$ and consider different $`N=(\mathrm{}_1,\mathrm{}_2,n)`$ and set $`F_\lambda ^sf_N^{\mu _0}`$ in equation (6). Further we omit $`\mu _0`$ from $`f_N^{\mu _0}`$ and write the relevant coupled set of equations (6) as

$$\left[\frac{d^2}{d\rho ^2}+1\frac{\nu _N(\nu _N+1)}{\rho ^2}\right]f_N^s+\underset{N^{}}{}\frac{2\alpha _{NN^{}}^s}{P\rho }f_N^{}^s=0.$$
(9)
For our numerical computations we truncate the set to some maximum value $`N_{mx}`$ of N. These $`N_{mx}`$ equations in $`N_{mx}`$ variables are needed to be solved from origin to infinity. Actually we need construction of $`N_{mx}`$ independent solutions which vanish at the origin. Now for convenience we divide the whole solution domain $`(0,\mathrm{})`$ into three subdomains $`(0,\mathrm{\Delta }),(\mathrm{\Delta },R_{\mathrm{}})`$ and $`(R_{\mathrm{}},\mathrm{})`$, where $`\mathrm{\Delta }`$ has the value of a few atomic units and $`R_{\mathrm{}}`$ is a point in the asymptotic domain. Best choices for these may be made by simple variations. Results do not depend significantly on these. But for converged results in some situations, values of $`R_{\mathrm{}}`$, as in H$``$M-SOW calculation , are to be thousands of atomic units. Next we proceed for solutions over subdomains. For $`(R_{\mathrm{}},\mathrm{})`$ we have simple analytic solutions :

$$f_{snN}^s(\rho )=\underset{\mathrm{}}{}\frac{a_{kN}^{(\mathrm{})}sin\theta _k}{\rho ^{\mathrm{}}}+\frac{b_{kN}^{(\mathrm{})}cos\theta _k}{\rho ^{\mathrm{}}}$$
(10)
$$f_{snN}^s(\rho )=\underset{\mathrm{}}{}\frac{c_{kN}^{(\mathrm{})}sin\theta _k}{\rho ^{\mathrm{}}}+\frac{d_{kN}^{(\mathrm{})}cos\theta _k}{\rho ^{\mathrm{}}}$$
(11)
where $`f_{snN}^{(k)}`$ and $`f_{csN}^{(k)}`$ are the N-th element of the k-th solution vectors. Obviously these give 2$`N_{mx}`$ independent solution vectors. The coefficients in these expressions are determined through recurrence relations (see Das ) in terms of $`a_{kN}^{(0)}=a_{kN}`$ and $`b_{kN}^{(0)}=0`$, $`c_{kN}^{(0)}=0`$, $`d_{kN}^{(0)}=a_{kN}`$, $`a_{kN}`$ being the N-th element of the k-th eigen vector of the charge matrix $`A=(\alpha _{NN^{}})`$. Here we have $`\theta _k=\rho +\alpha _kln\mathrm{\hspace{0.33em}2}\rho `$, $`\alpha _k`$ being the k-th eigen value of A.

Solution over $`(\mathrm{\Delta },R_{\mathrm{}})`$ is also very simple. Because of the simple structure of equations (9) a Taylors expansion method works nicely. Earlier for the (e, 2e) problem, Das also adopted this approach . But the main difficulty lies in the construction of the solution vectors over $`(0,\mathrm{\Delta })`$. In those calculations on (e, 2e) problem Das used an approach as in R-matrix calculations . But very often, this invites pseudo resonance type behaviour causing undesirable oscillations in the cross sections. Recently  for (e, 2e) problem, we applied the finite difference method (a five-point scheme) for solutions in the interval $`(0,\mathrm{\Delta }`$) and thereby get rid of undesirable oscillations. So we adopted here the same approach but with a seven-point scheme in place of the five-point scheme.

Thus for the solution in the interval $`(0,\mathrm{\Delta })`$ we recast equations (9) in terms of R instead of $`\rho `$, as

$$\left[\frac{d^2}{dR^2}+P^2\frac{\nu _N(\nu _N+1)}{R^2}\right]f_N^s+\underset{N^{}=1}{\overset{N_{mx}}{}}\frac{2\alpha _{NN^{}}^s}{R}f_N^{}^s=0,$$
(12)
and solve these equations as a two point boundary value problem by difference equation method. At $`R=0`$, the solution vectors are set to zero while at $`R=\mathrm{\Delta }`$ we assign to the k-th solution vector the k-th column of the unit matrix. The matrix for the corresponding difference equation is a sparse matrix and for its solution special methods are available. Here we use biconjugate gradient method . We find that this method readily works and gives converged solutions.

Now for the difference equations we divide the interval $`[0,\mathrm{\Delta }]`$ into m subintervals of length h with mesh points

$$0=R_0<R_1<R_2<\mathrm{}<R_k<\mathrm{}<R_{m1}<R_m=\mathrm{\Delta }$$

with $`R_k=R_0+kh`$ and use the following seven-point difference formula:

$`f_N^{^{\prime \prime }}(R_k)`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}[{\displaystyle \frac{1}{90}}f_N(R_{k3}){\displaystyle \frac{3}{20}}16f_N(R_{k2})`$ (13)
$`+`$ $`{\displaystyle \frac{3}{2}}f_N(R_{k1}){\displaystyle \frac{49}{18}}f_N(R_k)+{\displaystyle \frac{3}{2}}f_N(R_{k+1}){\displaystyle \frac{3}{20}}16f_N(R_{k+2})`$
$`+`$ $`{\displaystyle \frac{1}{90}}f_N(R_{k+3})]+\{{\displaystyle \frac{69}{25200}}h^6f_N^{(viii)}(\xi _1)\}`$
for $`k=3,4,\mathrm{},m4,m3`$, and for $`k=1,2`$ and $`m2,m1`$ the formulae

$`f_N^{^{\prime \prime }}(R_1)`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}[{\displaystyle \frac{3}{8}}f_N(R_0)+6f_N(R_1){\displaystyle \frac{11}{2}}h^2f_N^{^{\prime \prime }}(R_2)`$ (14)
$``$ $`{\displaystyle \frac{51}{4}}f_N(R_3)h^2f_N^{^{\prime \prime }}(R_3)+6f_N(R_4+{\displaystyle \frac{3}{8}}f_N(R_4)]`$
$`+`$ $`\{{\displaystyle \frac{23}{10080}}h^6f^{(viii)}(\xi _2)\}.`$
$`f_N^{^{\prime \prime }}(R_2)`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}[{\displaystyle \frac{3}{8}}f_N(R_1)+6f_N(R_2){\displaystyle \frac{11}{2}}h^2f_N^{^{\prime \prime }}(R_3){\displaystyle \frac{51}{4}}f_N(R_3)`$ (15)
$``$ $`h^2f_N^{^{\prime \prime }}(R_4)+6f_N(R_4)+{\displaystyle \frac{3}{8}}f_N(R_5)]`$
$`+`$ $`\{{\displaystyle \frac{23}{10080}}h^6f^{(viii)}(\xi _3)\}.`$
$`f_N^{^{\prime \prime }}(R_{m2})`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}[{\displaystyle \frac{3}{8}}f_N(R_{m5})+6f_N(R_{m4})h^2f_N^{^{\prime \prime }}(R_{m4}){\displaystyle \frac{51}{4}}f_N(R_{m3})`$ (16)
$`{\displaystyle \frac{11}{2}}h^2f_N^{^{\prime \prime }}(R_{m3})+6f_N(R_{m2})+{\displaystyle \frac{3}{8}}f_N(R_{m1})]`$
$`+\{{\displaystyle \frac{23}{10080}}h^6f^{(viii)}(\xi _4)\}.`$
and

$`f_N^{^{\prime \prime }}(R_{m1})`$ $`=`$ $`{\displaystyle \frac{1}{h^2}}[{\displaystyle \frac{3}{8}}f_N(R_{m4})+6f_N(R_{m3})h^2f_N^{^{\prime \prime }}(R_{m3}){\displaystyle \frac{51}{4}}f_N(R_{m2})`$ (17)
$`{\displaystyle \frac{11}{2}}h^2f_N^{^{\prime \prime }}(R_{m2})+6f_N(R_{m1})+{\displaystyle \frac{3}{8}}f_N(R_m)]`$
$`+\{{\displaystyle \frac{23}{10080}}h^6f^{(viii)}(\xi _5)\}.`$
The quantities on the right hand sides within curly brackets represent the error terms. The corresponding difference equations are obtained by substituting these expressions the values of second order derivatives from the differential equation (12). For continuing these solutions in the domain $`(\mathrm{\Delta },R_{\mathrm{}})`$ we need first order derivatives $`f_{}^{}{}_{N}{}^{}(R)`$ at $`\mathrm{\Delta }`$. These are computed from the difference formula

$`f_N^{^{}}(R_m)`$ $`=`$ $`{\displaystyle \frac{1}{84h}}[f_N(R_{m4}+24f_N(R_{m2})128f_N(R_{m1})+105f_N(R_m)]`$ (18)
$`+`$ $`{\displaystyle \frac{2h}{7}}f_N^{^{\prime \prime }}(R_n)+\{{\displaystyle \frac{4h^4}{105}}f_N^{(v)}(\xi )\}`$
Here too, the quantity within curly brackets represents the error term. The solutions thus obtained in $`(0,\mathrm{\Delta })`$ are then continued over $`(\mathrm{\Delta },R_{\mathrm{}})`$ by Taylor’s expansion method, as stated earlier, with stabilization after suitable steps . The $`N_{mx}`$ independent solution vectors so obtained, are put together to get the solution matrix $`f_0`$. The solution matrices $`f_{sn}`$ and $`f_{cs}`$ are similarly obtained, whose N-kth element are respectively $`f_{snN}^{(k)}`$ and $`f_{csN}^{(k)}`$, given by (10) and (11) respectively.

Next we introduce the K-matrix through the relation

$$f_0B=f_{sn}+f_{cs}K$$
(19)
where $`B`$ is an unknown constant matrix. The K-matrix is determined from matching values and first order derivatives at $`R_{\mathrm{}}`$, where all of $`f_0,f_{sn}`$ and $`f_{cs}`$ are valid. (It may be noted here that there is a slight departure in our definition of K-matrix from the usual practice. However, it is symmetric as it should be).

Finally the physical scattering state with appropriate boundary conditions is taken as

$$𝐟_{\mathrm{𝐩𝐡}}=f_0\text{g}$$
(20)
and also we have

$`𝐟_{\mathrm{𝐩𝐡}}`$ $`=`$ $`(f_{sn}+f_{cs}K)𝐜`$ (21)
$`=`$ $`f_{sn}𝐜+f_{cs}𝐝`$
with

$$𝐝=K𝐜.$$
(22)
Thus the physical state is completely determined once the vector $`𝐜`$ is determined. Now $`𝐜`$ is determined from the consideration that $`\mathrm{\Psi }_{fs}^{()}`$ is asymptotically a (distorted) plane wave (representing the two outgoing electrons) plus incoming waves only. So the coefficients of the outgoing wave $`exp(i\rho )`$ of both $`\mathrm{\Psi }_{fs}^{()}`$ and the symmetrized plane wave (equation (5)) must be the same (except for the distorting term $`exp(i\alpha _kln\mathrm{\hspace{0.33em}2}\rho )`$). This requires

$$𝐜=[I+iK]^1𝐏$$
(23)
where

$$𝐏=2e^{i\frac{\pi }{4}}X^1𝚽^𝐬(\omega _\mathrm{𝟎}),$$
(24)
and X is the matrix comprising of the columns of eigen vectors of the charge matrix A and $`𝚽^𝐬`$ is given by

$$𝚽^𝐬(\omega _\mathrm{𝟎})=\left(\begin{array}{c}\varphi _1^s(\omega _0)\\ \mathrm{}\\ \varphi _{N_{mx}}^s(\omega _0)\end{array}\right).$$
(25)
Finally the DPI triple differential cross section is given by

$$\frac{d^3\sigma }{d\mathrm{\Omega }_1d\mathrm{\Omega }_2dE_1}=\frac{2\pi ^2\alpha p_1p_2}{\omega _i}|T_{fi}|^2,$$
(26)
after the inclusion of $`\mu _0`$-part of $`\mathrm{\Psi }_{fs}^{()}`$ in $`T_{fi}`$.

## 3 Results

In our present calculation we have applied the above hyperspherical partial wave approach both in length and velocity gauges. Here we consider equal-energy-sharing case, since then the computational problem becomes little simple. For this calculation we have chosen $`\mathrm{\Delta }=5`$ a.u., $`R_{\mathrm{}}=300`$ a.u., $`h=0.05`$ a.u. upto $`\mathrm{\Delta }`$ and 0.1 a.u. beyond $`\mathrm{\Delta }`$. We have included 90 coupled channels with $`n`$ upto 9 and $`(l_1,l_2)`$ combinations nearly as in ECS calculation  for electron \- hydrogen ionization collision. We have chosen the case of ionization at 20 eV excess energy as it has been widely considered and for which there are interesting experimental results . For the present calculations with 90 channels and $`R_{\mathrm{}}=300`$ a.u., our single differential cross section (SDCS) is little above the desired value of about 0.93 Kb/eV at $`E/2`$ (E being the excess photon energy). So we normalized our TDCS by scaling with a factor 0.8 (which is also the factor we use to scale our SDCS to get the desired value of 0.93 Kb/eV at $`E/2`$) both in the length and in the velocity gauges. The TDCS results thus obtained are presented in figure 1. Here we compare our results with the experimental results of Bräuning et al  and with the theoretical results of the CCC calculation only , since the overall agreement of the CCC results are known to be somewhat better compared to the results of other calculations like TDCC, H$``$M-SOW etc. In all the cases the agreement between the velocity and length gauge calculations is excellent everywhere. We also did calculations in the acceleration gauge and these are indistinguishable from those of the velocity gauge. Agreement with the experimental results and with the CCC results are also generally good except for some spurious peaks at $`\theta _{p_1}=0^o`$. The CCC results appear a little better compared to ours.

For unequal energy sharing our approach also works and we have reasonably good results. But for fully converged results, we have to consider larger values for $`R_{\mathrm{}}`$. So we wish to report such results in future when our study is completed.

## 4 Conclusions

The present calculation, reported here, has only approximately converged. The results we have obtained, go to show that the hyperspherical three-particle scattering state wave function, used in the present calculation, must be reasonably accurate from small distances to the asymptotic region, since the results in all the three gauges are practically identical. In contrast, the 3C or other similar wave functions, which are not accurate at finite distances, show strong gauge dependence . We also mention that the present calculation is free from any genuine difficulty and does not show any weakness worth mentioning. The present approach may easily be applied to DPI with varied types of polarization of the incident photons. With judicial choice of the parameter $`R_{\mathrm{}}`$ and possibly with the availability of better computational facilities, the method may be applied from very low energy to high energy cases. At this point if we recall the capability of the hyperspherical partial wave approach in representing electron-hydrogen-atom ionization collisions  at low energies (and also consider situations of very low energy cases, with excess energy 1 eV and below for which we had to take $`R_{\mathrm{}}`$ about 4000 - 5000 a.u. and get reliable results for ionization cross sections ), consider the present success, then we may expect the hyperspherical partial wave theory to have a very good prospect.

We are grateful to H. Bräuning for providing us with the experimental results and to Igor Bray and Anatoly Kheifets for providing us with the CCC results in electronic form. KC acknowledges support from the UGC in the form of a Minor Research Project F.PSW-035/02(ERO). SP is grateful to CSIR for providing a research fellowship.

## References

Figure Captions

Figure 1. Triple differential cross sections for photo double ionization of the helium atom for equal energy sharing geometry for 20 eV excess energy and for a) $`\theta _{p_1}=0^o`$, b) $`\theta _{p_1}=30^o`$, c) $`\theta _{p_1}=60^o`$, d) $`\theta _{p_1}=90^o`$, $`\theta _{p_1}`$ being measured from the photon polarization direction. Theory : continuous curve, present calculation in velocity gauge; dashed curve present calculation in length gauge; dotted curve, CCC calculation ; Experiment : absolute measured values of Bräuning et al.


text:
# Domino tilings and related models: space of configurations of domains with holes

## 1 Introduction

In the last ten years, a lot of progress has been done about the study of tilings. Most remarkably, W. P. Thurston , using work of J. H. Conway and J. F. Lagarias , introduced the notion of height functions, which encode domino tilings and calisson tilings of a polygon $`P`$.

The notion of height function appears to be a very powerful tool for the study of tilings. It has notably been extended by different authors   to study tiling algorithms for other sets of prototiles.

For domino tilings, height functions induce a lattice structure on the set of tilings of a fixed polygon (see ). Some important results are obtained from this structure: A linear time tiling algorithm , rapidly mixing Markov chains for random sampling  , computation of the number of necessary flips (local transformations involving two dominoes) to pass from a fixed tiling to another fixed tiling , efficient exhaustive generation of tilings  .

Dominoes are of particular importance to theoretical physicists, for whom dominoes are models of *dimers*, which are diatomic molecules (such as dihydrogen), and each tiling is seen as a possible state of a solid or a fluid.

The present paper tries to generalize previous results to figures which are not polygons, i. e. figures with holes. This is done by the introduction of an equilibrium function on edges of cells of the figure. With this tool, we prove that the set of tilings of any finite figure has a distributive lattice structure, of which we give a geometrical interpretation of this structure. To this end, we also need to introduce some structural notions, following works of J. Propp  and J. C. Fournier : The critical cycles, which induce forced components and generalized flips.

Our approach is constructive; this allows us to exhibit algorithms to compute the objects introduced. As a consequence, we obtain an exhaustive generation algorithm and a uniform random sampling algorithm.

We finish by proving that these ideas can be directly adapted for other types of tilings: Calisson tilings and tilings with bicolored Wang tiles.

## 2 Figures in the plane grid

### 2.1 The plane grid

Let $`\mathrm{\Lambda }`$ be the plane grid of the Euclidean plane $`^2`$. A vertex of $`\mathrm{\Lambda }`$ is a point whose coordinates are both integers.

A vertex $`v=(x_1,y_1)`$ is a *neighbour* of another vertex $`v=(x_2,y_2)`$ if $`|x_1x_2|+|y_1y_2|=1`$. Hence, each vertex $`v`$ has four neighbours $`v+(1,0)`$, $`v(1,0)`$, $`v+(0,1)`$ and $`v(0,1)`$ which are canonically called the East, West, North and South neighbour of $`v`$, respectively. An *edge* of $`\mathrm{\Lambda }`$ is the closed segment of straight line between two adjacent vertices. A cell of $`\mathrm{\Lambda }`$ is a (closed) unit square whose corners are vertices. Two cells are 4-neighbours (respectively 8-neighbours) if they share an edge (respectively at least a vertex).

A directed graph $`G=(V,E)`$ is *symmetric* if $`(v,v^{})E`$ if and only if $`(v^{},v)E`$ for all $`v,v^{}V`$. In this work we deal with the symmetric directed graph (denoted by $`\mathrm{\Lambda }^+`$) obtained from the planar grid $`\mathrm{\Lambda }`$ by replacing each edge $`vv^{}`$ by two *arcs* $`(v,v^{})`$ and $`(v^{},v)`$. For an arc $`a=(v,v^{})`$ of $`\mathrm{\Lambda }^+`$ we denote by $`[a]`$ its associated edge in $`\mathrm{\Lambda }`$.

A *(directed) path* $`P`$ in a directed graph $`G=(V,E)`$ is a sequence of vertices $`(v_0,\mathrm{},v_k)`$ such that $`(v_i,v_{i+1})`$ is an arc of $`G`$ for every $`i=0,\mathrm{},k1`$. We denote by $`E(P)`$ the multiset of all the arcs used by the path $`P`$ and by $`V(P)`$ the multiset of its vertices. We say that $`G`$ is *connected* if any two vertices of $`V`$ are linked by a path.

A path $`P=(v_0,\mathrm{},v_k)`$ with $`v_k=v_0`$ is called a cycle. The cycle is elementary if $`v_i=v_j`$ and $`ij`$ imply $`\{i,j\}=\{0,k\}`$. In a plane graph one has two kinds of elementary cycles: The *clockwise* cycles and the *counterclockwise* ones.

Let $`G=(V,E)`$ be a symmetric directed graph. A function $`g:E`$ is *skew-symmetric* if $`g(v,v^{})=g(v^{},v)`$ for every $`(v,v^{})E`$. Given any function $`h:V`$ we define its associated difference function $`D(h):E`$ by $`D(h)(v,v^{})=h(v^{})h(v)`$, for all $`(v,v^{})E`$. Conversely, if $`G`$ is connected, given a function $`g:E`$ which satisfies $`g(C)=0`$ for all cycle $`C`$ of $`G_F`$ and a vertex $`w_0`$ of $`V`$, there exists a unique function $`h:V`$ such that $`h(w_0)=0`$ and $`D(h)(a)=g(a)`$ for all $`aE`$.

Let $`E^{}`$ be a multiset of arcs of $`G`$. We denote by $`g(E^{})`$ the sum of the values $`g(a)`$ over all the arcs $`aE^{}`$ (each arc $`a`$ is counted according to its corresponding multiplicity in the multiset). Then $`g(E^{})=_{aE^{}}g(a)`$. For a path $`P`$, instead of $`g(E(P))`$ we use the shorthand $`g(P)`$.

We assume that cells of $`\mathrm{\Lambda }`$ are colored as a checkerboard. We thus have black cells and white cells, and two cells sharing an edge have different colors. Let us define the *spin* function $`sp`$ on the arcs of $`\mathrm{\Lambda }^+`$. For each arc $`a=(v,v^{})`$, the spin of $`a`$ is noted $`sp(a)`$ and given by:

* $`sp(a)=1`$ if an ant moving from $`v`$ to $`v^{}`$ has a white cell on its left side (and a black cell on its right side);
* $`sp(a)=1`$ otherwise.

For each clockwise elementary cycle $`C`$, one has $`sp(C)=4Dis(C)`$ where $`Dis(C)`$, the *disequilibrium* of $`C`$, is the difference between the number of black cells and the number of white cells enclosed by $`C`$. The result is true for each cycle enclosing a single cell, and each elementary cycle can be decomposed into such sqaure cycles.

### 2.2 Figures

A figure $`F`$ of $`\mathrm{\Lambda }`$ is a 4-connected, finite union of cells of $`\mathrm{\Lambda }`$. The unique infinite 8-connected component of $`^2F`$ is denoted by $`H_{\mathrm{}}`$. The other ones are called the holes of $`F`$. The set of all edges in the boundary of $`F`$ is denoted by $`E_b(F)`$. The set of edges in $`H_{\mathrm{}}F`$ is called the *outer-boundary* of $`F`$ and denoted by $`E_{ob}(F)`$. Analogously, we denote by $`V_b(F)`$ the set of all the vertices on the boundary of $`F`$.

Because of the two types of connectivity for cells, we replace (until the end of the paper) each vertex $`v`$ of $`F`$ such that each edge issued from $`v`$ is on the boundary of $`F`$, by two vertices $`v_1`$ and $`v_2`$, each of them connected to exactly two neighbours of $`v`$ (see Figure 1).

A figure defines a symmetric directed graph $`G_F=(V_F,E_F)`$ such that $`V_F`$ is the set of corners of cells of $`F`$ (once duplication is done), and $`E_F`$ is the set of arcs $`a`$ such that $`[a]`$ is a side of a cell of $`F`$. From this point of view, the clockwise and counterclockwise contours of each hole are elementary cycles of $`G_F`$.

For each elementary clockwise cycle of $`F`$ (*i. e.* whose arcs are in $`E_F`$) we define $`Dis_F(C)`$ as the difference between the number of black cells of $`F`$ and the number of white cells of $`F`$ enclosed by $`C`$.

### 2.3 Equilibrium function

Informally, we can say that we want to work as if $`F`$ had no hole. To this end, the informal idea is to introduce values on edges which make holes disappear. Precisely, this is done by the use of equilibrium functions as defined below:

###### Definition 2.1

An equilibrium function (denoted by $`eq`$) is a skew-symmetric function from $`E_F`$ to $``$ such that $`sp(C)+eq(C)=4Dis_F(C)`$, for every clockwise cycle $`C`$ of $`F`$.

For figures without holes, it suffices to take $`eq=0`$.

Notice that, through the decomposition of cycles, a function $`eq`$ is an equilibrium function if and only if the following conditions holds.

* $`eq(C)=0`$ for each elementary cycle $`C`$ around a cell of $`F`$.
* $`eq(C)=sp(C)`$ for each cycle $`C`$ which follows clockwise the boundary of a hole of $`F`$.

More generally, for every cycle $`C`$ of $`F`$, $`eq(C)`$ is a number which does not depend on the chosen equilibrium value. We prove in Section 4 that every figure has an equilibrium function which can be efficiently computed.

We also need some auxiliary functions deduced from the function $`eq`$.

###### Definition 2.2

The functions $`eq_r`$, $`𝐭`$ and $`𝐛`$ are defined as follows.

1. $`eq_r(a)=eq(a)sp(a)`$ for all $`aE_F`$.
2. $`𝐭(a)=eq(a)sp(a)+2`$ for all $`aE_FE_b(F)`$.
3. $`𝐛(a)=eq(a)sp(a)2`$ for all $`aE_FE_b(F)`$.
4. $`𝐭(a)=𝐛(a)=eq(a)+sp(a)`$ for all $`aE_b(F)`$.

Note that for any arc $`a`$ in $`E_F`$, $`𝐭(a)𝐛(a)`$ is either 0 or 4.

## 3 The lattice of tilings

In this section we associate three classes of objects to a figure: Tilings, height functions and acyclic orientations. Our goal is the study of tilings, and height functions and acyclic orientations are some powerful tools to do this study.

### 3.1 Tilings and height functions

A domino is a figure formed by two cells sharing an edge, which is called the central axis of the domino. A tiling $`T`$ of a figure $`F`$ is a set of dominoes included in $`F`$, with pairwise disjoint interiors (*i. e.* there is no overlap), such that the union of the tiles of $`T`$ equals $`F`$ (*i. e.* there is no gap). Each tiling $`T`$ of a figure $`F`$ is completely determined by the set of its central axis. The characteristic function of a tiling $`T`$ defined from the set of arcs of $`F`$ into $``$ is given by $`\chi _T(a)=1`$ if $`[a]`$ is a central axis of a tile of $`T`$ and $`\chi _T(a)=0`$ otherwise.

Let $`T`$ be a tiling of $`F`$. The height difference in $`T`$, noted $`g_T`$ , is the skew-symmetric function defined by

$$aE_Fg_T(a)=eq_r(a)+2sp(a)(12\chi _T(a))$$

Then $`g_T(a)\{𝐛(a),𝐭(a)\}`$ for every $`aE_F`$.

Let us define $`g_F(a)`$ as $`eq(a)+sp(a)`$ for every $`aE_F`$. It can be seen that for each pair $`(T,T^{})`$ of tilings of $`F`$ and every $`aE_F`$, $`g_T(a)g_T^{}(a)=4sp(a)(\chi _T(a)\chi _T^{}(a))`$. Thus, if $`g_T=g_T^{}`$, then $`T=T^{}`$. Hence, the function $`g_T`$ is a tool to encode the tilings. Moreover, for every arc $`aE_b(F)`$ and every tiling $`T`$ of $`F`$ we necessarily have $`g_T(a)=g_F(a)`$. Thus $`g_T(a)`$ does not depend on $`T`$ for $`aE_b(F)`$. Additionally, $`g_T(a)g_T^{}(a)\{4,0,4\}`$, for every $`aE_F`$.

###### Proposition 3.1

Let $`T`$ be a tiling of a figure $`F`$. For each cycle $`C`$ of $`F`$, one has $`g_T(C)=0`$.

This proposition is a generalization of a theorem by J. H. Conway  about tilings of polygons.

* Proof. (sketch) It suffices to prove the result for elementary cycles since the height difference of each cycle is the sum of the height differences of the elementary cycles which compose it. This is done by induction on the number of cells of $`\mathrm{\Lambda }`$ enclosed by the cycle.

The case of a cycle following the boundary of a hole is easily treated from the definition of equilibrium functions. We also verify that the proposition holds for elementary cycles of length 4 around a cell.

We now use induction. If we are not in one of the cases treated above, then the area enclosed by the cycle can be cut by a path in $`F`$, which induces two new cycles, each of them enclosing less cells of $`\mathrm{\Lambda }`$ than the original cycle. Thus, by the induction hypothesis, the height difference of both induced cycles is null, from which it is easily deduced that the height difference of the original cycle is null.

Proposition 3.1 guarantees the correctness of the definition below.

###### Definition 3.2

For each tiling $`T`$, the height function induced by $`T`$ (denoted by $`h_T`$) is the function from the set $`V_F`$ of vertices of cells of $`F`$ (once necessary vertex duplications have been done) to the set $``$ of integers, defined by $`h_T(w_0)=0`$ and $`D(h_T)=g_T`$.

We now give a characterization of height functions of tilings.

###### Proposition 3.3

Let $`w_0`$ be a fixed vertex of $`H_{\mathrm{}}`$. We denote by $`_F`$ the class of all the functions $`h:V_F`$ satisfying the following properties:

* $`h(w_0)=0`$;
* for each arc $`a`$ of $`E_F`$, $`D(h)(a)\{𝐛(a),𝐭(a)\}`$.

For each tiling $`T`$, the function $`h_T`$ belongs to $`_F`$. Conversely, for each $`h_F`$ there exists a tiling $`T`$ such that $`h=h_T`$.

* Proof. The first statement follows directly from the definition of $`g_T`$. Let $`h`$ belong to $`_T`$ and let $`C`$ be a cycle around a cell. Then $`D(h)(C)eq_r(C)=sp(C)`$ since $`D(h)(C)=eq(C)=0`$. Clearly, $`|sp(C)|=4`$. Since $`D(h)(a)\{𝐛(a),𝐭(a)\}`$, there are three arcs of $`C`$ such that $`D(h)(a)eq_r(a)=2sp(a)`$ and exactly one arc satisfying $`D(h)(a)eq_r(a)=2sp(a)`$. Thus, the set $`T`$ of all dominoes whose central axis $`[a]`$ is such that $`hd(a)eq_r(a)=2sp(a)`$ is a tiling $`T`$ of $`F`$. The equality $`h(v)=h_T(v)`$, for each vertex $`v`$ of $`F`$, is obvious by induction on the distance from $`w_0`$ to $`v`$.

The proposition above allows one to consider each tiling as a height function.

###### Lemma 3.4

For every pair of height functions $`h`$ and $`h^{}`$ and for each vertex $`v`$ of $`F`$, one has $`h(v)h^{}(v)=0[4]`$. Moreover, $`h(v)h^{}(v)`$ does not depend on the chosen equilibrium function.

* Proof. Obvious by induction on the length of a shortest path from $`w_0`$ to $`v`$.

Let $``$ be the canonical order on functions: $`hh^{}`$ if and only if $`h(v)h(v^{})`$ for all $`vV_F`$.

###### Proposition 3.5

Let $`h`$ and $`h^{}`$ belong to $`_F`$. The functions $`inf(h,h^{})`$ and $`sup(h,h^{})`$ belong to $`_F`$.

In the vocabulary of order theory (see for example , ) the above proposition can be restated as follows: $`(_F,)`$ is a distributive lattice.

* Proof. Let $`h_1=inf(h,h^{})`$. We shall prove that for every arc $`(v,v^{})E_F`$, $`h_1(v^{})h_1(v)\{h(v^{})h(v),h^{}(v^{})h^{}(v)\}`$ (the proof for $`sup(h,h^{})`$ is similar).

For the sake of contradiction, let us assume that there exists an arc $`a=(v,v^{})`$ of $`F`$ such that $`h_1(v)=h(v)<h^{}(v)`$ and $`h_1(v^{})=h^{}(v^{})<h(v^{})`$. From Lemma 3.4, one has $`h^{}(v^{})h(v^{})4`$ and $`h(v)h^{}(v)4`$. Then $`\alpha :=h^{}(v^{})h(v^{})+h(v)h^{}(v)`$ satisfies $`\alpha 8`$. On the other hand, $`\alpha =h^{}(v^{})h^{}(v)(h(v^{})h(v))`$. Since $`h`$ and $`h^{}`$ belong to $`_F`$, one obtains $`\alpha 𝐛(v,v^{})𝐭(v,v^{})=4`$ which contradicts the hypothesis.

We define the following order on $`\mathrm{\Gamma }_F`$: $`TT^{}`$ if $`h_Th_T^{}`$. From Proposition 3.3, $`(\mathrm{\Gamma }_F,)`$ is isomorphic to $`(_F,)`$. Thus $`(\mathrm{\Gamma }_F,)`$ is a distributive lattice.

### 3.2 Forced components

###### Definition 3.6

An elementary cycle $`C`$ of $`G_F`$ is critical if $`𝐭(C)=0`$.

Such a cycle is strongly critical if, moreover, $`sp(a)=1`$, for each arc $`aE(C)E_b(F)`$.

We say that $`v`$ and $`v^{}`$ are *critically equivalent* if there is exists a critical cycle $`C`$ such that $`v,v^{}V(C)`$. The equivalence classes of this equivalence relation are called the *forced components* of the figure $`F`$.

Let $`\widehat{G}_F=(\widehat{V}_F,\widehat{E}_F)`$ be the symmetric graph whose vertices are the forced components and where $`(U,U^{})`$ is an arc of $`\widehat{G}_F`$ if there exists $`vU`$ and $`v^{}U^{}`$ such that $`(v,v^{})`$ is an arc of $`G_F`$. This graph is called the *graph of forced components* of $`F`$.

Notice that each boundary cycle of a hole of $`F`$ is strongly critical. The boundary cycle of the outer boundary is critical if and only if $`F`$ contains as many black cells as white cells. A strongly critical cycle can be deduced from each critical cycle by replacing each interior arc $`a`$ with $`sp(a)=1`$ by a sequence $`a^{},a^{\prime \prime },a^{\prime \prime \prime }`$ of three arcs of positive spins.

It is easy to see that for any cycle $`C`$ of $`G_F`$, $`𝐭(C)`$ and $`eq(C)`$ do not depend on the chosen equilibrium function. Thus, the notion of critical cycle only depends on the shape of the figure.

###### Definition 3.7

Let $`T`$ be a tiling of $`F`$. The graph of $`T`$ is the spanning subgraph of $`G_F`$, denoted by $`G_T=(V_F,E^{})`$, where $`aE^{}`$ if and only if $`g_T(a)=𝐭(a)`$.

By definition, every arc in $`E_b(F)`$ is an arc of $`G_T`$. Moreover, $`aE_TE_b(F)`$ if and only if $`\chi _T(a)=0`$ and $`sp(a)=1`$, or $`\chi _T(a)=1`$ and $`sp(a)=1`$.

The following proposition is the reason why we are interested in critical cycles and forced components.

###### Proposition 3.8

Let $`C`$ be a cycle of $`G_F`$.

* If $`C`$ is critical, then for every tiling $`T`$ the cycle $`C`$ is a cycle of $`G_T`$. Conversely, if $`C`$ is a cycle of $`G_T`$ for some tiling $`T`$, then $`C`$ is critical.
* If $`C`$ is strongly critical, then for every tiling $`T`$ the cycle $`C`$ is a cycle of $`G_T`$ which does not cut any tile of $`T`$.

* Proof. If $`C`$ satisfies $`𝐭(C)=0`$, then for every tiling $`T`$ of $`F`$ one has $`𝐭(C)=g_T(C)`$. Therefore $`𝐭(a)=g_T(a)`$ for every $`aE(C)`$, whence $`C`$ is a cycle of $`G_T`$. Conversely, if $`C`$ is a cycle of $`G_T`$, then $`g_T(a)=𝐭(a)`$ for every $`aE(C)`$. Thus, from Proposition 3.1, $`𝐭(C)=g_T(C)=0`$ and $`C`$ is a critical cycle.

For the second part, if $`C`$ is strongly critical then $`𝐭(a)>g_T(a)`$ implies $`sp(a)=1`$, from which one knows that $`\chi _T(a)=0`$; this means that $`C`$ does not cut any tile of $`T`$. Conversely, let $`[a]`$ be an interior edge which does not cut any tile of $`T`$. By definition, $`g_T(a)=sp(a)+eq(a)`$. Moreover, since, $`a`$ belongs to $`E_T`$, one has $`g_T(a)=𝐭(a)=eq(a)sp(a)+2`$ and finally $`sp(a)=1`$.

###### Corollary 3.9

If $`F`$ has a strongly critical cycle $`(v_0,v_1,\mathrm{},v_p)`$ such that for each integer $`0i<p`$, $`[v_i,v_{i+1}]`$ an interior edge, then there exists no tiling of $`F`$.

* Proof. Let $`v_j=(x_j,y_j)`$ be the vertex of this cycle with $`x_j+y_j`$ maximal, and, moreover, $`x_j`$ minimal with respect to the previous condition. One necessarily has $`v_{j1}=v_j+(1,0)`$, $`v_{j2}=v_{j1}+(0,1)`$ and (moreover) $`v_{j+1}=v_j+(0,1)`$. Now, follow the cycle until a vertex $`v_{j+2k}`$ such that $`v_{j+2k}v_j+(k,k)`$ (see Figure 3).

Let $`T`$ be a tiling of $`F`$; at least one tile of $`T`$ must be cut by an edge of the path formed from the part of the cycle from $`v_j`$ to $`v_{j+2k}`$. But this is impossible, from Proposition 3.8. Thus there exists no tiling.

From Corollary 3.9, if $`F`$ can be tiled then there are three kinds of forced components: The component $`U_{\mathrm{}}`$, the single components which are reduced to a single vertex, and the hole components which contain the contour of at least one hole.

The following lemma establish an useful relation between the height functions and the forced components.

###### Lemma 3.10

Let $`v`$ and $`v^{}`$ be critically equivalent vertices. For all $`h`$ and $`h^{}`$ in $`_F`$, one has $`h(v)h^{}(v)=h(v^{})h^{}(v^{})`$.

* Proof. Let $`C=(v_0,v_1,\mathrm{},v_p)`$ be a critical cycle passing through $`v`$ and $`v^{}`$. One can assume without loss of generality that $`v=v_0`$ and $`v^{}=v_k`$. One has $`h(v^{})h(v)=_{i=0}^{k1}D(h)(v_i),(v_{i+1}))=_{i=0}^{k1}𝐭(v_i,v_{i+1})=h^{}(v^{})h^{}(v)`$, which yields the result.

Let us choose one vertex $`v_U`$ in each forced component $`U`$ of $`G_F`$. From Lemma 3.10, $`hh^{}`$ if and only if $`h(v_U)h^{}(v_U)`$ for all $`U\widehat{V}_F`$.

We define a distance on $`_F`$ by

$$\mathrm{\Delta }(h,h^{}):=\underset{U\widehat{V}_F}{}\left|h(v_U)h^{}(v_U)\right|$$

Notice that the distance satisfies the following equalities:

$$\begin{array}{c}\mathrm{\Delta }(h,h^{})=\mathrm{\Delta }(h,inf(h,h^{}))+\mathrm{\Delta }(inf(h,h^{}),h^{})\hfill \\ \mathrm{\Delta }(h,h^{})=\mathrm{\Delta }(h,sup(h,h^{}))+\mathrm{\Delta }(sup(h,h^{}),h^{})\hfill \end{array}$$

### 3.3 Acyclic orientations and flips

In this part we prove that the lattice of height functions of a tileable figure $`F`$ is isomorphic to a lattice of a subclass of orientations of the graph of forced components. These lattices have been precisely studied by J. Propp .

###### Definition 3.11

A directed graph $`G`$ is an *orientation of $`F`$* if $`G`$ is an orientation of $`\widehat{G}_F`$ such that $`|E(C)E(G)|=\frac{1}{4}𝐛(C)`$ for all cycle $`C`$ of $`\widehat{G}_F`$. We denote by $`𝒢_F`$ the class of all the acyclic orientations of $`F`$.

Let $`h`$ be a height function and let $`T`$ be the corresponding tiling. Let $`G_h=(\widehat{V}_F,E)`$ be the *graph of strongly connected components* of $`G_T`$, *i. e.* the directed graph defined by $`(U,U^{})E`$ if and only if there exists an arc of $`G_T`$ from a vertex of $`U`$ to a vertex of $`U^{}`$.

It is well known (see  for example) that the graph of strongly components of any directed graph is acyclic. Hence, $`G_h`$ is an acyclic orientation of $`\widehat{G}_F`$.

###### Proposition 3.12

Let $`h`$ be an element of $``$. Then $`G_h`$ belongs to $`𝒢_F`$. Conversely, for each $`G=(V,E)𝒢_F`$ there exists $`h`$ such that $`G=G_h`$.

* Proof. Let $`C`$ be a cycle of $`\widehat{G}_F`$. Then $`0=D(h)(C)=D(h)(E(C)E)+D(h)(E(C)E)`$. Since $`D(h)(E(C)E)=𝐭(E(C)E)`$ and, for any arc $`a`$ of $`E_F`$, $`𝐭(a)=𝐛(a)+4`$, one obtains: $`D(h)(E(C)E)=𝐛(E(C)E)+4|E(C)E|`$. Moreover, from the definition of $`G_h`$ one obtains $`D(h)(E(C)E)=𝐛(E(C)E)`$. Finally, $`0=4|E(C)E|+𝐛(C)`$.

Conversely, let $`G=(V,E)`$ belong to $`𝒢_F`$. Let $`g`$ be the function defined by $`g(a)=𝐭(a)`$ if $`aE`$ and $`g(a)=𝐛(a)`$ if $`a\widehat{E}_FE`$. We prove that $`g(C)=0`$ for all cycle $`C`$. Clearly $`g(C)=g(E(C)E)+g(E(C)E)`$. By definition of $`g`$, one has $`g(E(C)E)=𝐭(E(C)E)`$ and $`g(E(C)E)=𝐛(E(C)E)`$. Since for any arc $`a`$ of $`E_F`$ $`𝐭(a)=𝐛(a)+4`$, one obtains that $`g(E(C)E)=𝐛(E(C)E)+4|E(C)E|`$. Then $`g(C)=𝐛(E(C)E)+4|E(C)E|+𝐛(E(C)E)=𝐛(C)+4|E(C)E|=0`$. Thus there exists $`h`$ such that $`h(w_0)=0`$ and $`D(h)=g`$.

###### Definition 3.13

Let $`G=(V,E)`$ belong to $`𝒢_F`$ and let $`UU_{\mathrm{}}`$ be in $`V`$ without incoming (resp. outgoing) arcs. The graph obtained from $`G`$ by an *upward* (resp. *downward*) *flip in $`U`$* is the acyclic directed graph $`G^U=(V,E^+)`$ (resp. $`G^U=(V,E^{})`$ ) where

$$E^+=E\{(U,U^{}):(U,U^{})E\}\{(U^{},U):(U,U^{})E\}$$

and

$$E^{}=E\{(U^{},U):(U^{},U)E\}\{(U,U^{}):(U^{},U)E\}$$

An upward or downward flip in $`U`$ corresponds to reversing all the arcs incident to $`U`$. If $`U`$ is reduced to a single vertex $`v`$, the flip is said local. If $`U`$ contains the contour of a hole, we say that it is a hole flip.

###### Proposition 3.14

For any $`G=(V,E)𝒢_F`$, $`G^U=(V,E^{})`$ belongs to $`𝒢_F`$.

* Proof. Let $`C`$ be a cycle of $`\widehat{G}_F`$. Since the number of arcs $`(U,U^{})`$ used by $`C`$ is equal to the number of arcs $`(U^{},U)`$ used by $`C`$, one has $`|E(C)E|=|E(C)E^{}|`$.

Let us denote by $`h^U`$ the unique function in $`_F`$ such that $`G_{h^U}=G^U`$. We have $`h^U=h`$ in $`VU`$ and $`|h^Uh|=4`$ in $`U`$. The two corresponding tilings differ only around $`U`$.

In particular, in $`\mathrm{\Gamma }_F`$ a local flip in $`U=\{v\}`$ (see Figure 4) is the replacement in $`T`$ of the pair of dominoes which cover the $`2\times 2`$ square centered in $`v`$ by the other pair which can cover the same square.

The upward flips defined above canonically induce an order on the set $`𝒢_F`$. Given $`G`$ and $`G^{}`$ in $`𝒢_F`$, we say that $`G_{flip^+}G^{}`$ if and only if there exists a sequence $`(U_0,\mathrm{},U_{p1})`$ of vertices of $`\widehat{G}_F`$ and a sequence $`(G_0,G_1,\mathrm{},G_p)`$ of graphs of $`𝒢_F`$ such that $`G_0=G`$, $`G_p=G^{}`$ and, for each integer $`0i<p1`$, $`G_{i+1}`$ is deduced form $`G_i`$ by an upward flip.

###### Proposition 3.15

Let $`h`$ and $`h^{}`$ belong to $`_F`$. Then $`h^{}h^{}`$ if and only if $`G_h_{flip^+}G_h^{}`$. Moreover, in this case, one can pass from $`G_h`$ to $`G_h^{}`$ by a sequence of $`\mathrm{\Delta }(h,h^{})/4`$ flips.

* Proof. The direct part of the proposition is proved by induction on the quantity $`\mathrm{\Delta }(h,h^{})`$. The result is obvious if $`\mathrm{\Delta }(h,h^{})=0`$ (i. e. $`h=h^{}`$).

Now, let us assume that $`\mathrm{\Delta }(h,h)0`$ and $`hh^{}`$. We will prove that there exists a forced component $`U`$ such that, for each vertex $`v_U`$ of $`U`$, $`h(v_U)<h^{}(v_U)`$ (which implies $`h(v_U)h^{}(v_U)4`$) and an upward flip can be done from $`G_h`$ on $`U`$.

Let $`U_0`$ be a component such that $`h(v_{U_0})<h^{}(v_{U_0})`$. If an upward flip can be done from $`G_h`$ on $`U_0`$, then we are done. Otherwise, there exists an arc $`(v_0,v_1)`$, with $`v_0`$ in $`U_0`$ and $`v_1`$ in another forced component $`U_1`$ such that $`D(h)(v_0,v_1)=𝐛(v_0,v_1)`$. From Lemma 3.10, we have: $`h^{}(v_{U_1})h(v_{U_1})=h^{}(v_1)h(v_1)`$ and $`h^{}(v_{U_0})h(v_{U_0})=h^{}(v_0)h(v_0)`$. Moreover, $`h^{}(v_1)h(v_1)=D(h^{})(v_0,v_1)+h^{}(v_0)h(v_0)D(h)(v_0,v_1)`$. Thus $`h^{}(v_{U_1})h(v_{U_1})=h^{}(v_{U_0})h(v_{U_0})+D(h^{})(v_0,v_1)𝐛(v_0,v_1)`$, which yields $`h^{}(v_{U_1})h(v_{U_1})h^{}(v_{U_0})h(v_{U_0})>0`$.

Either an upward flip can be performed on $`U_1`$, or the same argument can be repeated from $`U_1`$ to obtain an arc $`(v_1,v_2)`$ from $`U_1`$ to another forced component $`U_2`$, $`D(h)(v_1,v_2)=𝐛(v_1,v_2)`$. By repeating the process, there are two possibilities: Either one obtains a forced component $`U_i`$ on which no upward flip can be done, or an infinite sequence $`(U_i)_i`$ is obtained. But since $`F`$ is finite, the second possibility would imply that there exists a finite subsequence $`(U_j,U_{j+1},\mathrm{},U_k)`$ such that $`U_j=U_k`$, which is a contradiction since $`G_h`$ is acyclic.

We have now proved the existence of a forced component $`U`$ such that an upward flip can be done from $`G_h`$ around $`U`$ to obtain a function $`h^U`$. Notice that $`\mathrm{\Delta }(h^U,h^{})=\mathrm{\Delta }(h,h^{})4`$. This proves that by induction one can pass from $`h`$ to $`h^{}`$ with $`\mathrm{\Delta }(h,h^{})/4`$ flips.

The direct part of the first part of the proposition is obvious.

###### Corollary 3.16

The function $`\phi :_F𝒢_F`$ defined by $`\phi (h)=G_h`$ is an order isomorphism between $`(_F,)`$ and $`(𝒢_F,_{flip^+})`$.

###### Corollary 3.17

Let $`h`$ and $`h^{}`$ belong to $`_F`$. The number of successive flips needed to pass from $`G_h`$ to $`G_h^{}`$ is $`\mathrm{\Delta }(h,h^{})/4`$. Moreover, the components $`U`$ on which flips are done in such a sequence are those such that $`h(v_U)h^{}(v_U)`$.

* Proof. A flip changes $`\mathrm{\Delta }(h,h^{})/4`$ by one unit. Thus $`\mathrm{\Delta }(h,h^{})/4`$ is a lower bound and the bound is reached if each flip lets the quantity decrease. Thus the bound can be reached only if the components $`U`$ on which flips are done are precisely those such that $`h(v_U)h^{}(v_U)`$.

Conversely, we have seen that for $`hh^{}`$, one can pass from $`G_h`$ to $`G_h^{}`$ by a sequence of $`\mathrm{\Delta }(h,h^{})/4`$ flips. For the general case we use $`inf(h,h^{})`$. Passing through $`G_{inf(h,h^{})}`$, one can pass from $`G_h`$ to $`G_h^{}`$ by a sequence of $`\mathrm{\Delta }(h,inf(h,h^{}))/4+\mathrm{\Delta }(inf(h,h^{}),h^{})/4=\mathrm{\Delta }(h,h^{})/4`$ successive flips.

###### Corollary 3.18

For each pair of tilings $`(T,T^{})`$, it is possible to pass from $`T`$ to $`T^{}`$ by a sequence of local flips if and only if $`h_T=h_T^{}`$ for all the vertices on the boundary of $`F`$. Moreover, in this case the number of local flips needed is $`_{vE_F}\left|h_T(v)h_T^{}(v)\right|`$.

* Proof. This is a special case of the previous corollary.

#### 3.3.1 Freeness and rigidity of tilings

Informally, our tools allow one to see what is forced and what one has to choose in order to tile a figure.

###### Lemma 3.19

Let $`h`$ belong to $`_F`$. The function $`h`$ is maximal if and only if for each forced component $`U`$ there exists a path $`P`$ from $`U_{\mathrm{}}`$ to $`U`$ in $`G_h`$.

The function $`h`$ is minimal if and only if for each forced component $`U`$ there exists a path $`P`$ from $`U`$ to $`U_{\mathrm{}}`$ in $`G_h`$.

* Proof. If there exists $`U^{}`$ such that there is no path from $`U_{\mathrm{}}`$ to $`U^{}`$ in $`G_h`$, then by taking a longest $`P^{}`$ which finishes in $`U^{}`$ one clearly deduces the existence of a component $`UU_{\mathrm{}}`$ of $`G_h`$ with no incoming arcs. Thus an upward flip can be done in $`U`$ and $`h`$ is not maximal.

Conversely, if $`hh_{\mathrm{max}}`$ then it is possible to pass from $`h`$ to $`h_{\mathrm{max}}`$ by a sequence of upward flips. The first component $`u`$ on which a flip is done in such a sequence has no outgoing arc.

The proof for $`h_{\mathrm{min}}`$ is similar.

###### Proposition 3.20

For each component $`U`$ such that $`UU_{\mathrm{}}`$, one has $`h_{\mathrm{min}}(v_U)h_{\mathrm{max}}(v_U)`$.

* Proof. Assume that $`h_{\mathrm{min}}(v_U)=h_{\mathrm{max}}(v_U)`$. Take an sequence of upward flips from $`h_{\mathrm{min}}`$ to $`h_{\mathrm{max}}`$. Since no flip is done on $`U`$, no upward flip can be done on each component on a path of $`G_{h_{\mathrm{min}}}`$ from $`U`$ to $`U_{\mathrm{}}`$. Thus this path is also a path of $`G_{h_{\mathrm{max}}}`$. By concatenating it with a path of $`G_{h_{\mathrm{max}}}`$ from $`U_{\mathrm{}}`$ to $`U`$, a cycle appears in $`G_{h_{\mathrm{max}}}`$, which is a contradiction.

###### Corollary 3.21

A vertex $`v`$ of $`V_F`$ belongs to $`U_{\mathrm{}}`$ if and only if $`h(v)=h^{}(v)`$ for all $`h,h^{}_F`$.

An arc $`a`$ links two vertices of the same component if and only if $`D(h)(a)=D(h^{})(a)`$, for all $`h,h^{}_F`$.

* Proof. If $`v`$ belongs to $`U_{\mathrm{}}`$, then $`h(v)h^{}(v)=h(w_0)h^{}(w_0)=0`$ from Lemma 3.10. Otherwise, $`h_{\mathrm{min}}(v)h_{\mathrm{max}}(v)`$ according to the previous Proposition.

If an arc $`a`$ links two vertices of the same component, then $`D(h)(a)=D(h^{})(a)`$ from Lemma 3.10. Conversely, if an arc $`a`$ does not link two vertices of the same component, then one of its vertices is in a component $`U`$ such that $`UU_{\mathrm{}}`$. Take an sequence of upwards flips from $`h_{\mathrm{min}}`$ to $`h_{\mathrm{max}}`$. A flip of the sequence is done on $`U`$, which implies that there exists a pair $`(h,h^{})`$ of functions such that $`D(h)(a)D(h^{})(a)`$.

## 4 Effective construction and algorithms

### 4.1 forced components

When $`F`$ can be tiled (which is the interesting case), the graph of forced components can be constructed in polynomial time: Given a tiling $`T`$, we have to construct the graph of strongly connected components of $`G_T`$, which can be done in linear time (see  or ).

We know from matching theory that there exists a $`O(n^{3/2})`$ algorithm  to obtain such a tiling $`T`$, where $`n`$ denotes the area of $`F`$. N. Thiant  gives an algorithm which is linear in the area enclosed by $`F`$ (i. e. the sum of the area of $`F`$ and the areas of the holes).

### 4.2 Construction of an equilibrium function

An equilibrium function can be exhibited using cut lines (see also ) as follows (see Figure 5):

For each hole $`H_i`$ of $`F`$, we (arbitrarily) fix a vertical segment $`L_i=[p_i,p_i^{}]`$ (which is called a cut line issued from $`H_i`$) of $`^2`$ such that $`p_i`$ is the central point of a highest cell of $`H_i`$; there exists a positive integer $`n_i`$ such that $`p_i^{}=p_i+(0,n_i)`$, the vertex $`p_i^{}`$ is not in $`F`$, and, for each integer such that $`0<n_i^{}<n_i`$, the point $`p_i+(0,n_i^{})`$ is the central point of a cell of $`F`$. Hence, the point $`p_i^{}`$ is the central point of a cell of another 8-connected component, $`H_j`$ of $`^2F`$, with $`ji`$ (and, possibly, $`j=\mathrm{}`$).

We say that $`H_j`$ is the (immediate) predecessor of $`H_i`$. This construction yields a directed tree whose vertices are 8-connected components of $`^2F`$. This tree is rooted in $`H_{\mathrm{}}`$. We inductively define the step value of a hole $`H_i`$ (denoted by $`step(i)`$) as follows: Let $`C_i=(v_{i,0},v_{i,1},\mathrm{},v_{i,p_i})`$ be a cycle which (clockwise) follows the boundary of $`H_i`$. We state:

$$step(i)=\underset{j|H_j\text{ has }H_j\text{ for predecessor}}{}step(j)sp(C_i)$$

We can now define the equilibrium value of an arc $`(v,v^{})`$ by:

* $`eq(v,v^{})=step(i)`$ if $`v^{}`$ is the East neighbour of $`v`$ and the line segment $`[v,v^{}]`$ crosses the cut line $`L_i`$;
* $`eq(v,v^{})=step(i)`$ if $`v^{}`$ is the West neighbour of $`v`$ and the line segment $`[v,v^{}]`$ crosses the cut line $`L_i`$;
* $`eq(v,v^{})=0`$ if the line segment $`[v,v^{}]`$ crosses no cut line.

This function satisfies the conditions of the definition: It is obvious for cycles surrounding single cells, and if $`C_i`$ is a cycle which clockwise surrounds the hole $`H_i`$, one has $`eq(C_i)=step(i)_{j|H_j\text{ has }H_j\text{ for predecessor}}step(j)=sp(C_i)`$.

The equilibrium function defined above has a specific property which can be used for algorithmic arguments: For each pair $`(v,v^{})`$ of vertices of $`V_F`$, there exists a path $`P`$ from $`v`$ to $`v^{}`$ such that for each arc $`a`$ of $`P`$, $`eq(a)=0`$. This yields that there exists a spanning tree $`T_F`$ of $`G_F`$ such that for each arc $`a`$ of $`T_F`$ (seen as a symmetric graph), $`eq(a)=0`$.

For each arc $`a`$ of $`E_F`$, $`|eq(a)|4n`$ (this upper bound can be reached by taking a cycle which uses $`a`$ and cuts no line, except $`a`$). A precise study shows that such an equilibrium function can be constructed in $`O(n\mathrm{log}n)`$ time units, where $`n`$ denotes the number of cells of $`F`$.

### 4.3 Minimal tiling

For figures without holes, W. P. Thurston  has exhibited an algorithm that builds the minimal tiling. For the general case, this algorithm can be generalized as follows.

Initialization: For each vertex $`v`$, the algorithm uses a variable value $`h(v)`$ and two fixed values $`low(v)`$ and $`sup(v)`$.

From previous results, we know that for all $`aE_{ob}(F)`$ and all tiling $`T`$, $`g_T(a)=g_F(a)`$. Then the value $`h(v)`$ can be computed for all the vertices on the boundary of $`H_{\mathrm{}}`$. If a contradiction appears, then stop. Otherwise, set $`low(v)=sup(v):=h(v)`$ for all $`v`$ on the boundary of $`H_{\mathrm{}}`$.

We first construct a spanning tree $`T_F`$ on $`G_F`$ rooted in $`w_0`$. Then, for all vertices $`v`$ not on the boundary of $`H_{\mathrm{}}`$, we set $`sup(v)=𝐭(T_Fv)`$ and $`low(v)=𝐛(T_Fv)`$, where $`T_Fv`$ denotes the unique path from $`w_0`$ to $`v`$ in $`T_F`$.

We set $`h(v)=𝐛(v)`$ for all vertices $`v`$ not on the boundary of $`H_{\mathrm{}}`$.

The algorithm also uses a set $`V`$ consisting of vertices $`v`$ such that there exists a neighbour $`v^{}`$ of $`v`$ satisfying $`h(v)+top(v,v^{})<h(v^{})`$. This set is computed during the initialization.

Main loop: While $`V`$ is not empty:

* Pick a vertex $`v`$ in $`V`$ and update $`h(v)`$ by adding 4 units. If $`h(v)>sup(v)`$ after updating, then stop (there is no tiling).
* Update $`V`$ by adding the neighbours $`v^{\prime \prime }`$ of $`v`$ such that $`h(v^{\prime \prime })+𝐭(v^{\prime \prime },v)<h(v)`$ and remove $`v`$ if necessary.

###### Proposition 4.1

Given a figure $`F`$ formed of $`n`$ cells, the above algorithm stops after at most $`n^2`$ time passages through the loop.

Moreover, the algorithm stops with $`V`$ empty if and only if there exists a tiling. In this case, when the algorithm stops, one has $`h=h_{T_{min}}`$.

* Proof. First remark that, for each arc $`(v,v^{})`$ of $`E_F`$, $`h(v^{})h(v)=𝐭(v,v^{})`$ in $`/4`$. This is true during the initialization: We introduce the cycle $`C`$ formed by the concatenation of (the opposite path of) $`T_Fv`$, $`T_Fv^{}`$ and $`(v,v^{})`$, we have $`𝐭(C)=0[4]`$ and this, together with the relations between $`𝐭`$ and $`𝐛`$, gives the result. Moreover, this property is preserved by the loop.

Thus each passage through the loop makes the sum $`_{vV_F}(sup(v)h(v))`$ decrease by at least 4 units since, on vertex $`v`$, $`sup(v)low(v)4|E(T_Fv)|4n`$ (where $`|E(T_Fv)|`$ denotes the number of arcs of $`T_Fv`$), so the algorithm stops after at most $`n^2`$ passages through the loop.

When $`V`$ becomes empty, each arc $`(v,v^{})`$ of $`E_F`$ satisfies the hypothesis of Proposition 3.3. Thus there exists a tiling $`T`$ such that $`h_T=h`$. Moreover $`hh_{T_{min}}`$ (this is true during the initialization, and this property is preserved by the loop), so $`h_{T_{min}}=h`$.

If the algorithm finds a vertex $`v`$ such that $`h(v)>sup(v)`$, then there is no tiling since, otherwise, for each tiling $`T`$, $`h_T(v)>sup(v)`$, which is a contradiction (clearly, from the definition of $`sup`$, $`h_T(v)sup(v)`$). This finishes the proof.

If we take into account the implementation, the algorithm’s cost is at most $`0(n^2)`$ time units, as follows:

All the values are encoded in unary numeration, which permits to add constant numbers in constant time. We use the equilibrium function described in part 4.2 and the spanning tree $`T_F`$ is chosen in such a way that for each arc $`a`$ of $`T_F`$, $`eq(a)=0`$. Thus, the initialization costs $`0(n^2)`$ time units, to compute $`h(v)`$ for each vertex and $`h(v^{})h(v)`$ for each arcs.

Each passage through the loop costs $`O(1)`$ time units since it consists in a fixed number of additions of 4 units and sign tests. This gives the time complexity.

Of course, a similar algorithm can be designed to construct the maximal tiling of $`F`$.

### 4.4 Exhaustive generation

An exhaustive generation can be done, extending ideas of  to figures with holes. Let $`(U_1,U_2,\mathrm{},U_q)`$ be a fixed total order of forced components of $`F`$ (except $`U_{\mathrm{}}`$). We define a total order $`<_{lex}`$ on tilings of $`F`$ as follows: Given two tilings $`T`$ and $`T^{}`$, we have $`T<_{lex}T^{}`$ if there exists an integer $`1ip`$ such that $`h_T(v)<h_T^{}(v)`$ for each vertex $`v`$ of $`U_i`$ and $`h_T(v)=h_T^{}(v)`$ for each vertex $`v`$ of $`U_j`$ with $`1j<i`$.

The order $`<_{lex}`$ is a linear extension of $`<_{height}`$, i. e. given two tilings $`T`$ and $`T^{}`$ such that $`T<_{height}T^{}`$, we have $`T<_{lex}T^{}`$.

###### Proposition 4.2

Let $`T`$ be a (non-maximal) tiling of $`F`$ and let $`T_{succ}`$ denote the successor of $`T`$ in the lexicographic order; let $`i`$ denote the largest integer such that an upward flip is possible in $`U_i`$.

The tiling $`T_{succ}`$ is the lowest tiling (for $`<_{height}`$) such that $`h_{T_{succ}}(v)=h_T(v)+4`$ for each vertex $`v`$ of $`U_i`$, and $`h_T(v)=h_{T_{succ}}(v)`$ for each vertex $`v`$ of $`U_j`$ with $`1j<i`$.

* Proof. Let $`T^{}`$ be a tiling such that $`h_{T_{succ}}(v)=h_T(v)+4`$ for each vertex $`v`$ of $`U_i`$ and $`h_T(v)=h_{T_{succ}}(v)`$ for each vertex $`v`$ of $`U_j`$ with $`1j<i`$. By definition, one has $`T_{succ}_{lex}T^{}`$.

Moreover, assume that $`h_T=h_{T_{succ}}`$ in $`U_i`$. Thus, by Corollary 3.16, one can pass from $`T`$ to $`sup(T,T_{succ})`$ by a sequence of upward flips, which contradicts the definition of the integer $`i`$.

This proposition enables one to generate all the tilings of $`F`$ as follows:

Initialization: Construct the graph $`G_{forc.}`$ of forced components, the tiling $`T_{min}`$, and $`o(T_{min})`$. Output the tiling $`T_{min}`$.

The algorithm uses a variable tiling $`T`$ stored in memory, which for initialization is equal to $`T_{min}`$.

Main loop: Compute the successor of $`T`$ as follows:

* Find the last component $`U_i`$ on which an upward flip can be done (if no upward flip is possible, then stop).
* Construct the minimal tiling $`T`$’ such that for each vertex $`v`$ of $`_{j=1}^{i1}U_j`$, $`h_T(v)=h_T^{}(v)`$ and for each vertex $`v`$ of $`U_i`$, $`h_T^{}(v)=h_T(v)+4`$.
* Replace $`T`$ by $`T^{}`$, $`o(T)`$ by $`o(T^{})`$, output the tiling and go back to the beginning of the loop.

The second item of the main loop can be done in $`O(n^2)`$ time units using an algorithm derived from the algorithm of construction of the minimal tiling (it suffices to change the initialization, fixing appropriate value of $`h(v)`$ for $`v`$ in $`_{j=1}^iU_j`$).

Thus, once the initialization is done, the maximal waiting time between two consecutive tilings is $`O(n^2)`$ time units. The memory space is $`O(n^2)`$ since for each vertex $`v`$, one has to store $`h_T(v)`$ (using unary numeration).

### 4.5 Uniform random sampling

Consider the following process: Given a tiling $`T`$, find at random a forced component $`C`$ and a direction (upwards or downwards). If a flip can be done in $`C`$ according to the chosen direction, then make this flip; otherwise, do not change $`T`$. Trivially, this Markovian random process is ergodic and converges to the uniform distribution.

Moreover, the method of “coupling from the past”  can be applied since the process is monotonic and one has a method to construct the maximal and minimal tilings. We thus have a randomized algorithm to sample domino tilings uniformly at random. The space required is polynomial.

It has been previously proved   that this process is rapidly mixing for figures without holes. We conjecture that is remains true in the general case.

About algorithms, the reader can also easily verify that, given a pair $`(T,T^{})`$ of tilings, one can

* compute in linear time if one can pass from $`T`$ to $`T^{}`$ by a sequence of local flips (it suffices to compare $`h_T`$ and $`h_T^{}`$ on the boundary of $`F`$);
* compute in polynomial time with (low degree) a shortest path of flips to go from $`T`$ to $`T^{}`$ (using the cyclic orientations and the distance) and the length of such a path.

## 5 Extension to other types of tilings

### 5.1 Calissons

The same study can easily be done for calisson (i. e. tiles formed by two neighbouring cells of the triangular lattice) tilings to get similar results. In this case, local flips are induced by the two tilings of hexagons formed by six triangular cells. There are only two small differences, detailed below:

* One has two types of connectivity for triangular cells (3-connectivity for cells which share an edge, and 12-connectivity for cells which share a vertex). Thus, some vertices have to be duplicated or triplicated (see Figure 6).
* For the proof of corollary 3.9, one has to consider a part of the critical cycle with vertices $`v=(x,y)`$ such that $`y`$ is maximal.

### 5.2 Bicolored Wang tiles

The case of dominoes is a particular case of tilings with Wang tiles (i. e. $`1\times 1`$ squares with colored edges, see  for details). They give rise to a tiling if the colors on the edges of neighbour squares are compatible.

An instance of the problem of tiling by bicolored Wang tiles is given by a finite figure and a coloration of the edges which are on its boundary. Hence, a domino tiling is a tiling by Wang tiles with one red edge and three blue edges, all the edges on the boundary of the figure being blue.

#### 5.2.1 Eulerian orientations

The same study can be done for tilings with “balanced Wang tiles” such that each square has two blue edges and two red edges. In this case, one has $`hd_T(v,v^{})=sp(v,v^{})+eq(v,v^{})`$ if the corresponding edge is blue and $`hd_T(v,v^{})=sp(v,v^{})+eq(v,v^{})`$ otherwise. This is recognizable as the height function for Eulerian orientations of the dual lattice, called the six-vertex ice model by physicists . This is also equivalent to the height function for three-colorings of vertices of the square lattice, and to alternating-sign matrices . The results are similar to those obtained for dominoes.

#### 5.2.2 Examples with finite height functions

For the case of “odd tiles” (i. e. tiles with exactly three edges of the same color (blue or red), see ): One has to take a height function in $`/8`$ such that $`hd_T(v,v^{})=sp(v,v^{})+eq(v,v^{})`$ if the corresponding edge is blue and $`hd_T(v,v^{})=3sp(v,v^{})+eq(v,v^{})`$ otherwise. With our technique of equilibrium value, it is easily proved that the set of the tilings of a fixed figure has a structure of boolean lattice (or hypercube), even if the figure has holes.

The case of “even tiles” (i. e. tiles with an even number of blue edges and an even number of red edges) is very similar, with $`val_T(v,v^{})=(sp(v,v^{}),0)`$ if the corresponding edge is blue and $`val_T(v,v^{})=(sp(v,v^{}),sp(v,v^{}))`$ otherwise. These values are taken in $`/4\times /2`$.