text:
# Dynamics and Topological Aspects of a Reconstructed Two-Dimensional Foam Time Series Using Potts Model on a Pinned Lattice

## 1 Introduction

This paper addresses the reconstruction of domain boundaries from experimental images of foams, cellular images which consist of sets of compact domains separated from each other by sharp boundaries which meet at vertices. The data to be reconstructed contain only partial information, in particular the positions of vertices (two dimensions) or edges (three dimensions) in the complete image. This task, which is a special case of the problem of image segmentation, recurs in many fields, e.g. medical imaging and cartography . An ISI Web of Science search reveals nearly 2000 references on image segmentation. Even when the image data is complete, a typical voxel image requires reconstruction to obtain information about domain shapes, volumes and adjacencies. In real images the presence of noise further complicates reconstruction .

The general reconstruction problem is underdetermined. Unless we supply further constraints, an arbitrary number of different domain patterns with different numbers of domains and domain area distributions can have the same set of vertices or edges. As a result, almost all reconstruction of three-dimensional foams has been manual, which is inconvenient and time consuming and greatly limits the numbers of bubbles which we can examine .

However, in the case of many domain patterns of real interest, e.g. metallic polycrystals, foams and biological cells, the patterns are well behaved in a manner which is clear to the eye, though difficult to quantify completely. For example, in two dimensions, boundaries meet at three-fold vertices at approximately 120 angles and are relatively smooth and constant in curvature. Many of these regularities, known as Plateau rules , result from the physical processes generating the domain boundaries: the presence of a surface tension or surface energy causes boundaries to be approximately minimal surfaces.

In these cases we can proceed quite far in reconstructing the complete two-dimensional pattern from the vertex positions and we present an algorithm that accomplishes this reconstruction with accuracy good enough for many practical purposes. The algorithm is attractive in that it uses a physical simulation of the development of minimal surfaces which in some sense duplicates the original formation of the pattern. One great advantage is that it works directly with real pixel-wise input data. We illustrate several pathological cases where the algorithm fails and discuss ways to fix these errors.

## 2 Motivation

Why do we need an algorithm that can reconstruct minimal surface partitions from vertex positions? Whole subfields of physics and mathematics have developed to study the structure of complex minimal surfaces (for a review see the web site of Kenneth Brakke and his program Surface Evolver at http://www.susqu.edu/facstaff/b/brakke/).

The division of space into subdomains with minimum partition area has attracted attention of scientists in fields ranging from Geography and Physics to Mathematics and Computer Science . The more restricted question of which geometric structure, replicated infinitely and filling space without leaving gaps, has minimum surface area, is known as the Kelvin Problem after Lord Kelvin who first proposed a solution in 1887 .

Dry liquid foams (i.e. foams in which the liquid fraction is below about 1%) are a good probe of the various proposed solutions to the Kelvin problem. The idealized model of a foam retains only the surface energy (proportional to the surface area) of the bubbles. Where bubbles contact each other, they form a cell face of negligible thickness. The bubble wall membranes bend, forming surfaces of constant mean curvature. The faces contact each other at edges known as Plateau borders where water accumulates and these edges in turn meet at Plateau vertices at the tetrahedral angle.

Foams also exhibit a complex dynamics in which gas passes from bubbles of higher pressure to lower pressure bubbles . Since smaller bubbles tend to have higher pressures because of geometric constraints, bubbles gradually disappear and the pattern coarsens . The same general phenomena with additional complications occur in metallic polycrystals, sintered ceramics, ferromagnets, micelles, ferrofluids and biological tissues and organs . Domain growth is a common feature of most of these patterns, though the time scales differ depending on the diffusion mechanisms controlling boundary motion. If the initial conditions are reasonably homogeneous, most of these materials develop geometrically similar structures . Thus, understanding the kinetics and geometry of a coarsening foam can provide information on a broad class of important materials.

Unfortunately, determining the three-dimensional structure of a foam is surprisingly difficult. Foams scatter light extremely effectively, so seeing deep inside a liquid foam is difficult . Optical tomography is only possible for very dry foams, in which the Plateau borders but not the faces are visible . Magnetic resonance imaging (MRI) can also reveal a foam’s three-dimensional structure . MRI measures the nuclear magnetic resonance (NMR) signal from protons in the hydrogen in the water present in the sample. However, since foams contain very little water, which mostly concentrates in the Plateau borders, the MRI signal is very weak and the best result we obtain is an image of the Plateau borders, which looks like a skeleton of the foam as in fig. 1. Recent synchrotron CT images provide higher quality images , but seeing the faces between bubbles remains effectively impossible.

The raw image is missing crucial information, e.g. the bubble volumes, because of the absence of the bubble membranes. The same problem can occur in confocal microscopy of cells, CT, MRI and ultrasound imaging of tissues and X-ray and neutron imaging of solids. An added difficulty is that the signal-to-noise ratio of such images is often close to one.

Reconstruction of foam structure from partial information is, then, essential to the study of cellular patterns. The concept underlying our reconstruction is energy minimization, in the sense that reconstructing the pattern is equivalent to finding the minimum energy partition consistent with the partial information available about the structure. Since the pattern itself results from constrained energy minimization, the use of an energy minimization method in reconstruction seems natural and appealing. In three dimensions, this partial information consists of a pixel image of the edges of bubbles (Plateau borders) and the reconstruction consists of finding the surfaces that connect the edges with minimum energy. In two dimensions, on the other hand, the corresponding partial data are the positions of the vertices of the domains and the reconstruction connects the vertices to obtain a minimum energy configuration. Although the results we present here concern two-dimensional reconstruction only, the method extends to three dimensions without major conceptual modifications.

We now summarize our method to reconstruct the two-dimensional structure starting solely from the positions of the vertices and our knowledge of minimal surfaces. Subsequent sections provide details and results. The method consists of three main independent steps:

* Finding the centers of bubbles in the pattern to be reconstructed.
* Reconstructing the pattern by nucleating a bubble from each center using the extended large-$`Q`$ Potts Model.
* Growing the nucleated bubbles until evolution stops.

Minimum distance maps provide the positions of the centers and the initial target areas (Each bubble receives a target area which we adjust during the simulations to achieve the minimum surface energy for each bubble) . A continuous pinning field traps the moving bubble edges.

We employ this method to reconstruct several patterns from a foam evolution time series obtained from a coarsening simulation. In order to test the validity of our results, we perform a detailed analysis of the dynamics and spatial structure of the reconstructed foam and compare it with the simulation and theoretical and experimental results found in the literature.

## 3 The Potts Model

The Potts Model  is a generalization of the Ising model to more than two spin components. Although initially proposed to study critical phenomena in statistical physics, the Potts Model finds a wide variety of applications, including the simulation of metallic grain growth , soap foams , magnetic froths  and biological cells .

We map the foam structure onto a rectangular lattice containing $`N_xN_y`$ sites. Each lattice site contains an integer $`\sigma `$ (a spin) which corresponds to a particular bubble (number of bubbles = number of spins). The boundary between two bubbles is the set of links between lattice sites associated to the spins of those bubbles. The boundary energy associates a positive energy (equivalent to a surface tension) with boundary links and zero energy for links within bubbles.

Coarsening has two time scales. Surface tension and viscous dissipation determine the time, $`\tau _1`$, which describes the relaxation of an edge towards equilibrium without volume (or area in two dimensions) variation. After a perturbation, the pattern returns to obeying the Plateau rules after a delay of $`\tau _1`$. The second time scale, $`\tau _2`$, depends mainly on gas permeability and foam polydispersity and controls the coarsening rate. In liquid foams $`\tau _1`$ is typically fractions of a second while $`\tau _2`$ is minutes or hours. In metals the two times are usually comparable, as they are in the simple Potts Model, so the foam only obeys Plateau’s rules approximately. Wejchert et al.  introduced an area constraint energy term to simulate equilibration in a foam, where boundary equilibration is much faster than diffusion. If we constrain the area for a fixed number of Monte Carlo steps (mcs as defined below) per volume diffusion step we can choose any ratio we like between the time scales.

The Hamiltonian including surface energy and the area constraint has the general form:

$$H=J\underset{i,j(sites)}{}(1\delta _{\sigma _i\sigma _j})+\lambda \underset{k(bubbles)}{}(S_kS_k^t)^2\text{.}$$
(1)
The first term in eq. 1 accounts for the surface energy. The sum in $`i`$ runs over all the sites in the lattice and the sum in $`j`$ runs over sites neighboring $`i`$. The parameter $`J`$ sets the energy/unit surface area. The second term is the area constraint. The sum in $`k`$ runs over all the bubbles in the pattern. $`\lambda `$ is a parameter specifying the strength of the area constraint, $`S_k`$ the current area of the $`k^{th}`$ bubble, and $`S_k^t`$ the target area of the same bubble. Because of the surface energy, each bubble’s area is usually smaller than its target area.

Most Potts simulation use a modified Metropolis Monte Carlo dynamics  in which we select a boundary lattice site at random and randomly propose to change its spin to the value of one of the neighboring spins (Kawasaki dynamics  and Glauber dynamics  are also possible). If the resulting change of energy is less than or equal to zero, we accept the new spin configuration. However, if the change in energy is positive, we accept the new configuration with Boltzmann probability:

$$P=\mathrm{exp}\left(\frac{\mathrm{\Delta }E}{kT}\right)\text{,}$$
(2)
where $`T`$ is the temperature assigned to allows for thermal fluctuations to overcome local energy minima. We define a Monte Carlo step as a sequence of $`N_xN_y`$ random site selections.

## 4 Generation of the input data

We may use any experimental picture of a two-dimensional foam as input data, once we know the vertex positions. To compare the statistical properties of the reconstructed foam to those of the original foam, we need a controlled series of well characterized foam structures. To develop and test our reconstruction procedures, we generate two-dimensional foam structures at different stages of a coarsening evolution using the Large-Q Potts Model .

The Hamiltonian contains the surface term only:

$$H=J\underset{i,j}{}(1\delta _{\sigma _i\sigma _j})\text{.}$$
(3)
The $`j`$ sum ranges up to third-nearest neighbors, or 20 sites around the $`i^{th}`$ site. In our simulation $`J=1`$. In order to minimize finite size and edge effects we use a $`1024\times 1024`$ square lattice with periodic boundary conditions. We start the simulations with $`16384`$ square $`8\times 8`$ bubbles and finish with just a few tens of bubbles.

We performed the main part of the numerical simulation at $`T=0`$. However, in order to overcome metastable traps, we periodically increased the temperature for short periods during the simulation, alternating 40 mcs at temperature $`T=0`$ (relaxation period) with 5 mcs at $`T=0.5`$ (fluctuation period). This method produces very realistic foam structures . Thus the Potts Model basis of both original and reconstructed structures does not bias our error estimates.

## 5 Pinning Field

Our inspiration for a continuous attractive pinning field to keep vertices fixed comes from Zener pins , which simulate the microstructural evolution of materials in the presence of a second phase dispersion of particles. The presence of insoluble precipitates inhibits grain growth, greatly affecting the mechanical properties of the material.

Srolovitz et al.  extended the Potts Model to simulate grain growth in the presence of a second phase particle dispersion. Their model incorporates particles by assigning a particular spin to the lattice sites the particles occupy. The kinetics is the regular Monte Carlo dynamics for the Potts Model. Other groups  have also used this procedure to obtain properties including grain size, characteristic exponents, etc. for such pinned materials. Krichevsky and Stavans  and Herrera et al.  performed two-dimensional experiments on foam evolution in the presence of pinning centers, using a square array of pins and randomly distributed pins respectively.

We construct a continuous attractive pinning field around a point size Zener pin at each vertex in the source image, contributing a Hamiltonian term:

$$H_p=\gamma \underset{i,j(sites)}{}(1\delta _{\sigma _i\sigma _j})\text{F}(x_i,y_i)=\gamma \underset{i,j(sites)}{}(1\delta _{\sigma _i\sigma _j})\times \left(\frac{1}{r_i}\right)\text{,}$$
(4)
where $`\gamma `$ is a constant that determines the pinning coupling strength and $`r_i`$ is the distance from the site $`i`$ to the closest pin. This term causes edges, and ultimately the reconstructed vertices, to move towards and eventually trap on the pins. The field F$`(x,y)`$ has its maximum value of $`1`$ at the position of the pins and falls off as the inverse of the distance ($`1/r`$). It contributes a negative energy at each place where a link between two different bubbles is present (see eq. 4). Therefore, boundaries and vertices cost less energy when located close to a pinning center, producing an attractive pinning force on the vertices.

We choose the value of $`\gamma `$ to produce a circular region around each center where the attraction of that center dominates, i.e., the center will definitely attract any edge or vertex within this circle. In our simulations, we chose the circles’ radii to be 20 pixels, which gives a value for the field, at the edge of the circles, of about 0.05. With $`\gamma =1000`$, the pinning energies are strong enough to ensure that almost no edges or vertices will slip off a pin. In the region between two pins closer than 40 lattice sites, the closest pin will attract the vertex due to the stronger pinning field in that direction. This situation becomes less significant when we impose the cutoffs discussed below.

Figure 2a shows the entire range of the pinning field for a particular foam configuration while fig. 2b shows the same field scaled to show the field within 20 pixels of the pins. The pinning field within the circles (visible in the top plane in fig. 2b) around each pin is much stronger than other energy components. The great advantage of this method compared to simple Zener pins is that in the final patterns almost all pins lie in edges. If the nucleation algorithm has missed bubbles, then the number of pins may be larger than the number of vertices in the reconstructed pattern and a few pins may not correspond to vertices. If the nucleation algorithm creates extra bubbles then a few vertices may not correspond to pins.

A drawback of this implementation is that the lattice becomes very stiff and, although the vertices pin appropriately, the edges are not free to relax. To overcome this problem, we impose a cutoff on the field, after a certain number of mcs. After the cutoff, the field still falls off as $`1/r`$ but we limit its range to a certain radius around each pin. This procedure keeps the vertices pinned but frees most of the edge length from the pinning field, allowing edges to straighten in order to minimize surface energy. We apply successively shorter range cutoffs until the range of the pinning field is about four lattice sites. Figure 2c-d shows the field configuration after a cutoff for $`r=`$ 20.

## 6 Reconstruction

### 6.1 Centers of Bubbles

Our first step is to find the bubbles’ centers. We follow Glantz and Prause and use a Euclidean distance map to find the centers . We construct such a map by picking a lattice site, calculating its distance from all the vertices and storing the minimum of these distances for that site. Repeating this procedure for all lattice sites produces the map . Figure 3a shows the map for a typical foam configuration. We assume that the centers of the bubbles sit on the local maxima in the map. Figure 3b shows the positions of the calculated centers (green) compared to the real ones (black) obtained from center of mass calculations. The agreement is fairly good. In three dimensions, this process is essentially the same, with the sole change that the distances are calculated from the sites to the edges.

This method, however, has some pathologies. Multiple maxima may be present within a bubble. We treat this problem by combining any two maxima that are closer to each other than to the closest vertex. The solution has an intrinsic flaw: when a bubble is long and thin, its centers may be closer to neighboring bubble centers than to the closest vertex, causing the disappearance of the bubble. Alternatively, the centers in a stretched bubble may be further from each other than from the nearest edge, causing the center map to create an extra bubble. Figure 4 illustrates these two situations.

Saddle points are another problem. The center point of an edge is a saddle point in the Euclidean map and if the edge is long enough, the saddle point may look like a maximum because of errors resulting from lattice discretization. An effort to exclude saddle points when searching for maxima may result in missing centers for small bubbles. If the initial bubble size distribution is narrow (bubbles are roughly the same size) and the domains are fairly isotropic, then no errors result. We are currently working to quantify the consequences and thresholds of these errors. Overall, we lose a few very small or elongated bubbles. During foam coarsening, these bubbles have a short lifetime (they disappear rapidly), so the error in most cases is quite small and has the effect of reconstructing a foam pattern that seems slightly later in time than the source pattern.

These pathologies result from the distance algorithm, i.e., the first part of our method. Our main reconstruction procedure, however, is independent of the method we use to find the centers of bubbles which we can improve or replace, e.g. the method proposed by Kammerer and Glantz .

### 6.2 Nucleation

Once we have chosen a set of bubble centers we nucleate a bubble with unique spin and area of one pixel at each center and set the remaining sites to spin 0. The spin 0 domain behaves like a big bubble which shrinks while the others coarsen. The bubbles coarsen under a target area constraint. These target areas are roughly the areas each bubble will have at the end of coarsening. We defined the pinning field above. The Hamiltonian is:

$$H=\underset{i,j(sites)}{}\left(J\frac{\gamma }{r_i}\right)(1\delta _{\sigma _i\sigma _j})+\underset{k(bubbles)}{}\mathrm{\Lambda }_k(S_kS_k^t)^2\text{,}$$
(5)
where we have generalized the area constraint. The area coupling constant $`\mathrm{\Lambda }_k`$ is:

$$\mathrm{\Lambda }_k=\lambda \left(\frac{1}{S_k^t}\right)S_{max}^t\text{,}$$
(6)
where $`\lambda `$ is a constant and $`S_{max}^t`$ is the largest of the target areas. With a uniform $`\mathrm{\Lambda }`$ as in previous applications, the area term is larger for bigger bubbles than for smaller bubbles. Therefore, the former ”swallow” the latter in order to minimize energy. To overcome this problem, we weight the coupling constant by making it inversely proportional to the bubble target area.

To estimate the tentative target area, we used the area of a circle centered at the bubble’s calculated center with radius equal to the distance to the closest vertex. Figure 3b illustrates this procedure. Clearly this approximation works better for isotropic bubbles. We could also use the average distance to the neighboring vertices or other estimates. The geometrical element that fills two-dimensional space with minimum surface energy is the hexagon. The area of a hexagon that fits inside a circle is about 10% smaller than the area of the circle itself. Therefore, we use 90% of the areas of the circles in fig. 3b as the tentative target areas. The method we use to obtain the tentative target areas plays only a small role, if any, in the final reconstruction.

To ensure that all the bubbles grow at comparable rates, we increase their target areas in steps to their full values. During the first 50 mcs we set the target areas to 10% of their full values, so after 50 mcs, all the bubbles have grown to 10% of their final sizes. After each 50 mcs, we increase the target areas by 10%, so the target areas reach their full values after 500 mcs.

After all the bubbles have achieved their final areas and the structure has pinned, we turn off the area constraint and impose a cutoff on the pinning field, as discussed above, to allow the interfaces to relax, and the areas to change, in order to obtain the lowest energy configuration.

For many purposes the approximate reconstruction is adequate. However, if we wish to measure the growth rates of bubbles to high accuracy we are examining small differences of large areas so even a small error in the areas gives a large error in their growth rates. We are developing procedures to correct the bubbles’ areas .

### 6.3 Results

Figure 5 shows (up to 25000 mcs) the growth of bubbles and structures reconstructed for a simulated foam coarsened for 50000 mcs. The last picture in the series shows a snapshot of the simulation superimposed on the reconstruction. Our reconstructions used the same alternating temperature approach as our coarsening simulations. We set $`T=1`$ to allow the vertices to find the pinning centers faster. The parameters for all reconstruction simulations were $`J=`$ 1, $`T=`$ 1, $`\gamma =`$ 1000, and $`\lambda =`$ 0.1. We repeated the same procedure for many different configurations and obtained the reconstructed foam coarsening of fig. 6.

### 6.4 Analysis of Results

During nucleation and growth of bubbles, the bubbles are at first circular and grow at the same pace due to the rescaling of their target areas. Eventually they impinge on one another, forming interfaces and vertices. The contact between bubbles keeps them from growing further and both the area constraint and pins further inhibit global coarsening.

Once the bubbles fill the lattice, their configuration is already close to the final one. Most of the vertices have already trapped on pinning centers and most untrapped vertices eventually trap during the relaxation of edges and rearrangement of areas after we turn off the target area constraint and impose cutoffs on the pinning field. The residual pinning field hinders further coarsening.

Although defects in the center calculation cause most errors in the reconstructed pattern, the nucleation of bubbles and relaxation of edges can introduce additional imperfections. A particular vertex may attach to the wrong pin or an edge may be missing. A bubble which was originally very long and thin might receive a target area much different from its original area because the anisotropy makes the nearest vertex much closer to the center than the average distance to the vertices. Since all bubbles nucleate as circles, an edge may miss a pin during the nucleation of a long, thin bubble, leaving a pin inside the bubble which never captures the edge. After this bubble achieves its final target area, it will still have to jump over an energy barrier in order to squeeze into a long, thin shape. A pair of stretched bubbles just before a $`T1`$ process (topological transition where a four-fold vertex decays into two three-fold vertices) in the original image will almost never reconstruct properly because it is so near the much lower energy state after the $`T1`$ that the reconstruction jumps over the real configuration to the lower energy configuration. However, the real foam pattern would evolve rapidly to this state so the effective error is small.

Even after a vertex pins, it may break loose. Vertices in less energetically favorable configurations, like those undergoing $`T1`$s or joining edges that make angles much smaller than 120, may detach from the pins despite very strong coupling to the pinning field. These extreme situations create loose vertices that evolve locally with the dynamics of an unpinned foam which may complete $`T1`$s in progress or lead to the disappearance of bubbles. As a result, the reconstructed structure may depart slightly from the original configuration as in fig. 6. Errors are less common in more homogeneous, isotropic foams. Figure 7 shows the fraction of missing bubbles as a function of time.

How do these deviations affect the global coarsening described in ? If we characterize coarsening by fitting the average area $`<A>`$ per bubble as a function of time to a power law:

$$<A>=kt^\alpha \text{,}$$
(7)
we obtain a growth exponent characterizing the coarsening. For a foam-like material, the characteristic exponent must be $`\alpha =1`$ . Besides the mean bubble area, the basic measures of the state of a foam are the distribution $`\rho (n)`$ of the number of sides (the probability that a randomly selected bubble has $`n`$ sides), and the normalized area distribution $`\rho (A/<A>)`$ (the probability that a bubble has an area which is a given fraction of the mean bubble area). The distributions show dynamical scaling, i.e. they become time invariant at long times. We define the $`m`$th moment of the side distribution as:

$$\mu _m=\underset{n=2}{\overset{\mathrm{}}{}}\rho (n)(n<n>)^m\text{,}$$
(8)
where $`<n>`$ is the average number of sides of a bubble in the pattern (for infinite patterns $`<n>=`$ 6). The second moment, $`\mu _2`$, measures the r.m.s. width of the distribution while $`\mu _3`$ and $`\mu _4`$ measure the asymmetry and the flatness of the side distribution.

To test these properties and quantify the errors, we analyzed the coarsening dynamics of both the simulated and reconstructed foams. Figure 8 shows the time evolution of $`<A>`$ for both foams. Both curves follow a power law, with exponents $`\alpha =1.03\pm 0.03`$ for the original foam (solid line in the figure) and $`\alpha =1.01\pm 0.03`$ for the reconstructed foam (dotted line). Both exponents are consistent, within error, with the expected value of one. Thus, the reconstructed foam has the same global dynamics as the original and theory.

Figure 9 shows the mean deviation of the reconstructed areas from the original areas given by:

$$<\frac{\mathrm{\Delta }A}{A}>=\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{1}{A_i^O}|A_i^OA_i^R|\text{,}$$
(9)
as a function of time, including and excluding missing bubbles in the calculation. The superscripts O and R stand for the original and reconstruction, respectively. The sum in eq. 9 runs over all $`N`$ bubbles in the foam. Both curves decrease slightly due to the increasing length scale of the pattern. The time averaged deviation is $`(9.2\pm 3.0)\%`$ including missing bubbles and $`(5.4\pm 1.6)\%`$ excluding them. Hence, the area of a reconstructed bubble differs by about 5.4% from the area of the original bubble.

This residual error is a consequence of reconstructing straight edges rather than the curved boundaries of a real foam. According to Plateau’s rules, the ratio between the perimeter and the square root of the area of a regular bubble is a constant ($`P/\sqrt{A}3.72`$) while the ratio for polygonal bubbles with straight sides is $`2\sqrt{n\mathrm{tan}(\pi /n)}`$. Therefore, the relative deviation from a regular Plateau bubble of a bubble with straight sides is approximately $`\delta (P/\sqrt{A})_n=2\sqrt{n\mathrm{tan}(\pi /n)}3.72`$. Figure 10 shows $`\delta (P/\sqrt{A})_n`$ (stars) on top of the time averaged relative difference between simulated and reconstructed areas given by:

$$<\frac{\mathrm{\Delta }A}{A}>_{time}=\frac{1}{N_n}\underset{i=1}{\overset{N_n}{}}\frac{1}{A_i^O}(A_i^OA_i^R)\text{,}$$
(10)
as a function of the number of sides $`n`$. The sum in eq. 10 runs over all $`N_n`$ bubbles with $`n`$ sides at any instant of time. For few-sided bubbles, the reconstructed area is smaller than the original area, while for many-sided bubbles the reconstructed area is larger than the original area. The reconstructed images had no 3- or 11-sided bubbles to analyze. Otherwise, the relative area deviations follow very closely the theoretical predictions for $`\delta (P/\sqrt{A})_n`$.

von Neumann’s law states that, in two dimensions, the area growth rate of individual bubbles with $`n`$ edges is proportional to $`(n6)`$:

$$\frac{dA_n}{dt}=\kappa (n6)\text{.}$$
(11)
Few-sided bubbles are convex, so a straight-sided reconstruction cuts off the bulges, while many-sided bubbles are concave so a straight-sided reconstruction fills in the dips. Figure 11 shows the deviation in von Neumann’s law calculated for the simulated and reconstructed foams as:

$$<\mathrm{\Delta }\left(\frac{dA_n}{dt}\right)>=\frac{1}{N_n}\underset{i=1}{\overset{N_n}{}}\left(\frac{dA_i^O}{dt}\right)^1\left|\frac{dA_i^O}{dt}\frac{dA_i^R}{dt}\right|\text{,}$$
(12)
where the summation in $`i`$ runs over all the $`N_n`$ $`n`$-sided bubbles at any instant of time. The nonlinearity seen in the inset to fig. 11 results from deviations in growth rate for small bubbles which are an artifact of the regular Potts Model simulation .

Figure 12 shows the time dependence of the areas of four individual bubbles. Filled and open symbols of the same type (circles, squares, etc.) correspond to simulated and reconstructed areas. The local dynamics are very similar for the simulation and reconstruction. Figure 6 shows that for the great majority of the bubbles, simulation and reconstruction follow the same dynamics.

Another important aspect of coarsening foam is self-similarity. Figure 13a shows the normalized area distributions at several times for the simulated foam. These distributions are congruent, indicating that the foam is in a scaling state. The thick line obtained by averaging all the distributions represents all the individual cases within $`3\%`$. The same applies to the reconstruction (fig. 13b). The scaled area distributions are equivalent within $`7\%`$ (fig. 14) so the reconstructed foam scales in the same way as the simulation.

Figure. 15 shows the evolution of the side distributions for both the simulation (a) and reconstruction (b). The distributions are time-independent and remarkably similar, peaking at $`n=6`$ . Figure 16 shows the second moment of the side distributions as function of time. The general behavior is similar for both patterns: $`\mu _2`$ oscillates about a stable value, since in the scaling state, the width of the distribution is nearly constant. The average values of $`1.62\pm 0.27`$ for the simulation and $`1.46\pm 0.30`$ for the reconstruction agree well with experimental values .

Figures 17 and 18 show, respectively, the third and fourth moments as functions of time. The third moment begins near its maximum ($`\mu _3=2.5`$) showing the inhomogeneity in the foam that lasts until about 30000 mcs. It then drops to a stable value of $`\mu _3=0.61\pm 0.27`$ for the simulation and $`\mu _3=0.63\pm 0.33`$ for the reconstruction. The fourth moment behaves similarly: an initial equilibration stage up to 30000 mcs, followed by stability with $`\mu _4=5.6\pm 1.4`$ for the simulation and $`\mu _4=5.5\pm 1.8`$ for the reconstruction, again agreeing with experiment .

## 7 Conclusions

Our method of reconstructing two-dimensional foam structures works well. The Potts Model generates minimal energy surfaces and the pins hinder coarsening. These two tools together reconstruct a coarsening foam. We show that imperfections mainly result from the distance map calculation, which is not the core of our method. A recent paper of Kammerer and Glantz suggests an alternative scheme which may eliminate many of these errors . The deviations do not induce systematic errors in the statistical properties of a foam during coarsening. The method extends to three dimensions with no major modifications, which will be very helpful in studies of the dynamics of three-dimensional cellular patterns.

Our aims are now two-fold. We must deal with real data. The reconstruction algorithm presented in this article will be the basis for analyzing the images of three-dimensional foams we recently obtained using X-ray tomography at the European Synchrotron Radiation Facility (ESRF). At this point, we stress that the major difference between two-dimensional and three-dimensional partial images is that the former consists of vertices while the later consists of edges of bubbles.

To extract the maximum information from experimental images we will directly couple the value of $`J`$ in the Potts Hamiltonian to the pixel gray level. A bubble boundary located in a dark region (fluid phase) will cost much less energy than a boundary in a light region (gas phase). This reformulation of the pinning field has the great advantage that it requires no image pretreatment or thresholding. Randomly distributed “fluid” pixels due to experimental noise will have much less power to pin boundaries than the organized pixels corresponding to edges, even if their gray level values are of the same order of magnitude. The method is also relatively immune to long wavelength inhomogeneities in the image saturation, e.g. due to MRI coil geometry.

We are working to improve our image reconstruction by developing heuristics to correct for missing/extra bubbles, by introducing additional mechanisms to deal with noise, and by implementing target area corrections.

We thank Dr. Mark Zajac and Dr. Marius Asipauskas for great help and infinite patience in our countless discussions, Dr. Igor Veretennikov and Dr. Burkhard Prause for the MRI image of a foam and Prof. François Graner for his helpful revisions to the manuscript. We acknowledge support from the Center for Applied Mathematics of the University of Notre Dame and NSF grants DMR-0089162 and INT98-02417, NASA grant NAG3-2366 and DOE grant DE-FGO299ER45785.


format: text

text:
# Lattice packings with gap defects are not completely saturated

## I Introduction

A packing of unit spheres in $`^n`$ is *$`k`$-saturated* if it is not possible to replace $`k1`$ spheres by $`k`$ and still have a packing; it is *completely saturated* if it is $`k`$-saturated for all $`k`$ FTK . As discussed in FTK , every completely saturated packing has maximum density, and every periodic packing with maximum density is completely saturated. Thus the densest lattice packing of circles in $`^2`$ (the *honeycomb* circle packing) is completely saturated. This motivates the question: Is the honeycomb circle packing the only one which is completely saturated?

G. Fejes Tóth, G. Kuperberg, and W. Kuperberg FTKK  asked whether a honeycomb circle packing with a linear shift defect, as in Figure 1, is completely saturated. This question appears to be closely related to the (still open) conjecture of L. Fejes Tóth LFT  asserting that the same packing is *solid*, meaning that if you obtain another packing by rearranging finitely many circles, it is congruent to the original one.

In this paper we will instead consider packings with a linear or planar gap defect in which the spheres on either side do not touch.

###### Definition.

Let $`𝒮`$ be a sphere packing in $`^n`$ let $`H^n`$ be a hyperplane, let $`H^+`$ be a closed half-space bounded by $`H`$, and let $`\stackrel{}{v}`$ be a vector perpendicular to $`H`$ in the direction of $`H^+`$. Let $`\stackrel{~}{𝒮}`$ be the packing obtained from $`𝒮`$ by moving all spheres with centers in $`H^+`$ by the vector $`\stackrel{}{v}`$. (Note that $`\stackrel{~}{𝒮}`$ is a packing because the motion does not decrease the distances between sphere centers.) Then $`\stackrel{~}{𝒮}`$ is a *sphere packing with a hyperplane gap defect* and $`\stackrel{}{v}`$ is the *width* of the gap.

Figure 2 shows a honeycomb circle packing with a linear gap defect.

###### Theorem 1.

A honeycomb circle packing in $`^2`$ with a linear gap defect is not completely saturated, regardless of the width and direction of the gap. An fcc (face-centered cubic lattice) packing of spheres in $`^3`$ with an arbitrary planar gap defect is also not completely saturated.

Since the honeycomb circle packing in $`^2`$ LA  and the fcc sphere packing in $`^3`$ GA  are the unique densest lattice sphere packings in 2 and 3 dimensions, we obtain the following corollary.

###### Corollary 2.

In $`^2`$ and $`^3`$, no lattice sphere packing with a hyperplane gap defect is completely saturated.

## II Honeycomb packing with gap defects

### II.1 Gap along a lattice line

Let $`𝒞`$ be the honeycomb circle packing in $`^2`$ with a circle center at every point of the form $`(2i+j,\sqrt{3}j)`$, with $`i,j`$. Let $`d>0`$ and move all circles with center on or above the $`x`$-axis by the vector $`\stackrel{}{v}=(0,d)`$ to obtain the packing $`\stackrel{~}{𝒞}`$.

We will show that $`\stackrel{~}{𝒞}`$ is not completely saturated. It suffices to rearrange finitely many circles in $`\stackrel{~}{𝒞}`$ to create a void with enough space for an extra circle. To achieve that, we will widen the gap by a factor bounded away from $`1`$ in a sufficiently large finite region. If we can create a gap of width $`d^{}=d+\delta `$, with $`\delta `$ a non-decreasing function of $`d`$, then by repeating the operation we can create a gap of any width. In particular, if the width is at least 2, there is room for another circle.

We define a *triangular block of size $`n`$* to be a set of $`n(n+1)/2`$ circles whose centers lie in an equilateral triangle of edge length $`2n2`$, *i.e.*, with $`n`$ circles along each edge. Note that a triangular block in $`\stackrel{~}{𝒞}`$ must lie entirely on one side of the gap. If the block points up, we will call it a $`\mathrm{\Delta }_n`$-block; otherwise we will call it a $`_n`$-block. The rearrangement to widen the gap is as follows: Take two $`\mathrm{\Delta }_n`$-blocks, two $`_n`$-blocks, and one $`\mathrm{\Delta }_{n1}`$-block that form a trapezoid, as in Figure 3. Move the two outer $`\mathrm{\Delta }_n`$-blocks down by a distance of $`d_1d`$, then move the two $`_n`$-blocks at $`60`$ degree angles from vertical by $`d_2`$, then move the middle $`\mathrm{\Delta }_n`$-block up by $`d_3`$, as in Figure 4. If we choose the distances $`d_1`$, $`d_2`$, and $`d_3`$ to maximize $`d_3`$, then $`d_3`$ is at least weakly monotonic in $`d`$, because increasing $`d`$ relaxes the constraints on the other parameters. Thus the gap is widened by $`\delta `$ monotonic in $`d`$, as desired. Since the parameter $`n`$ is arbitrary, it can be taken large enough to repeat the operation with smaller blocks to reach any desired width.

### II.2 Gap along an arbitrary line

Let $`𝒞`$ be the above honeycomb circle packing and let $`\mathrm{}`$ be any line, and now let $`\stackrel{~}{𝒞}`$ be the packing $`𝒞`$ with a gap along $`\mathrm{}`$. Our plan in this case is to create another (finite) gap parallel to a lattice line of $`𝒞`$ and reduce to the case of Section II.1. The gap can be made any length without sacrificing width. For each $`n3`$, we can choose a triangular block $`T`$ in $`𝒞`$ of length $`n`$ so that $`\mathrm{}`$ intersects the triangle formed by its circle centers, moreover so that $`\mathrm{}`$ does not meet the vertices. Then $`\mathrm{}`$ divides $`T`$ into two sub-blocks $`T_1`$ and $`T_2`$, one of which, say $`T_1`$, has only one corner of $`T`$. In $`\stackrel{~}{𝒞}`$, these two sub-blocks are separated by the gap, as in Figure 5. In $`\stackrel{~}{𝒞}`$, then, $`T_1`$ can be moved towards $`T_2`$ to introduce gaps on its other two sides. The length $`n`$ is arbitrary and the width of the new gaps does not depend on $`n`$. This completes the reduction.

## III The fcc packing with a gap

The fcc sphere packing in $`^3`$ decomposes into honeycomb layers, and, as noted by Kepler JK , it also decomposes into square-lattice layers. We prove below that the fcc-lattice sphere packing with a gap defect is never completely saturated, no matter how narrow the gap is and no matter along which plane the gap is formed. As before, we first establish special cases.

### III.1 Gap parallel to a square layer

Let $``$ be the fcc packing in $`^3`$ whose sphere centers are at $`(2i+k,2j+k,\sqrt{2}k)`$ with $`i,j,k`$. Let $`\stackrel{~}{}`$ be the packing $``$ with sphere centers in the upper half-space moved by $`\stackrel{}{v}=(0,0,d)`$ for some $`d>0`$. In this model the square-lattice layers are parallel to coordinate planes of $`^3`$, so $`\stackrel{~}{}`$ has a square-layer gap.

To show that $`\stackrel{~}{}`$ is not completely saturated, we will delete finitely many spheres, then rearrange finitely many others to create a void large enough to accommodate more spheres than the ones deleted. We define a *square pyramid block of size $`n`$* to be a collection of $`n(n+1)(2n+1)/6`$ spheres in $`\stackrel{~}{}`$ whose sphere centers lie in a square pyramid of edge length $`2n2`$, *i.e.*, with $`n`$ spheres along each of its edges. We define a *simplex block of size $`n`$* to be a collection of $`n(n+1)(n+2)/6`$ spheres in $`\stackrel{~}{}`$ whose spheres centers lie in a regular simplex of edge length $`2n2`$. Both kinds of blocks have $`n`$ spheres along each edge; see Figure 6 for examples.

By analogy with Section II.1, we will widen the gap by moving 5 square pyramid blocks and 4 simplex blocks forming a cross gable roof, as in Figure 7. The bases of the square pyramid blocks lie on the gap. We move the four outer square pyramid blocks down, then move the simplex blocks down and out, then move the middle square pyramid up after removing its top sphere. As before, we choose the motions in order to maximize the size of the last motion; in this arrangement the last motion is weakly monotonic in $`d`$. Thus the gap widens by $`\delta `$ with $`\delta `$ weakly monotonic in $`d`$. As before, $`\delta `$ does not depend on $`n`$.

As in Section II.1, $`n`$ can be large enough to iterate the procedure many times. Let $`k`$ be the number of iterations needed to widen the gap to 2; $`k`$ only depends on $`d`$. One sphere is deleted in each iteration, for a total of $`k`$ spheres removed. Since the void at the end grows with $`n`$, we can choose $`n`$ large enough to accommodate these $`k`$ spheres and then some.

### III.2 Gap parallel to a honeycomb layer

Our strategy in this case is to expose a square-layer gap in a large finite region by moving a large block that abuts the honeycomb gap. This reduces the problem to the case considered in Section III.1. The shape of the block is a cuboctahedron sawed in half to expose a regular hexagon face, as in Figure 8. As the figure indicates, the three square faces lie in square-lattice directions. Moving the block creates gaps along these faces.

### III.3 Gap parallel to an arbitrary plane

Assume now that the plane of the gap does not belong to any of the four parallel classes of the honeycomb layers. In a manner similar to that of Section II.2, we construct a large tetrahedron block which is cut by the gap into two sub-blocks $`T_1`$ and $`T_2`$. If the plane is not parallel to a square-lattice layer, then we can arrange that $`T_1`$ contains 1 vertex of the tetrahedron and $`T_2`$ contains 3 vertices. The sub-block $`T_1`$ can then be moved towards the opposite face, which creates a honeycomb-layer gap and reduces the problem to that of Section III.1.

## IV Open problems and conjectures

1. Is there a completely saturated sphere packing with a gap defect in $`^n`$ ($`n2`$)?
2. Let $`𝒮`$ be the densest lattice packing of spheres in $`^4`$, which is composed of fcc layers. Let $`H`$ be a hyperplane in an fcc direction and let $`\widehat{𝒮}`$ be a sphere packing obtained from $`𝒮`$ by introducing a gap defect along $`H`$. We conjecture that $`\widehat{𝒮}`$ is not completely saturated.
3. Let $`𝒯`$ be a tiling of $`^2`$ with regular hexagons, and let $`\widehat{𝒯}`$ be a packing obtained from $`𝒯`$ by introducing a gap. If the gap is sufficiently narrow, is $`\widehat{𝒯}`$ completely saturated?


text:
# CP symmetry and the strong interactions

## I Introduction

The SU(3) non-Abelian gauge theory of the strong interactions is quite remarkable in that, once an arbitrary overall scale is fixed, the only parameters are the quark masses. Using only a few pseudo-scalar meson masses to fix these parameters, the non-Abelian gauge theory describing quark confining dynamics is unique. It has been known for some time dashen  that, as these parameters are varied from their physical values, exotic phenomena can occur, including spontaneous breakdown of CP symmetry.

The possibility of a spontaneous CP violation is most easily demonstrated in terms of an effective chiral Lagrangian. In Section II I will review this model for the strong interactions with three quarks, namely the up, down, and strange quarks. This lays the groundwork for the discussion in Section III of the CP violating phase. Section IV discusses how heavier states, most particularly the $`\eta ^{}`$ meson, enter without qualitatively changing the structure.

Included among the parameters of the strong interactions is a complex phase which, if present, explicitly violates CP symmetry. This parameter appears to be extremely small dipole  since no such violation is seen phenomenologically. A puzzle for grand unification asks why is CP violation small for the strong interactions but not the weak pecceiquinn . It is sometimes suggested that a massless up quark would solve this problem, and I turn to this issue in section V. There I argue that the conventional view is incorrect. The concept of a single massless quark is renormalization scheme dependent and not physically meaningful. Indeed, due to non-perturbative effects, even the sign of a small up-quark mass is not uniquely defined. These effects are of higher order in the chiral expansion. But as long as the other quark masses are not identically zero, some ambiguity remains in the definition of the up-quark mass. This is likely not relevant for most phenomenological issues, but is unacceptable for solving something fundamental like the strong CP problem.

Because renormalization is required, the concept of an “underlying basic Lagrangian” does not exist in quantum field theory. Instead there are underlying symmetries, and the continuum theory is defined in terms of those and a few renormalized parameters. A single massless quark is not represented by any symmetry, and, because of confinement, its mass is not appropriate for a renormalized parameter. To elucidate these points, I review the non-perturbative definition of the continuum theory and the corresponding ambiguities in the quark mass. These issues remain even with the recently discovered chirally symmetric lattice fermions. Finally, Section VI contains some concluding remarks, including possible impacts of the CP violating structures for lattice gauge simulations.

## II The effective model

A CP violating phase appears naturally in the simplest chiral sigma model of interacting pseudo-scalar mesons. In this section I review the basic model and the standard connections between the quark masses and the meson masses. Nothing in this section is new; I am setting the stage for later discussion.

To be specific, consider the three flavor theory with its approximate SU(3) symmetry. Using three flavors simplifies the discussion, although the CP violating phase can also be demonstrated for the two flavor theory following the discussion in mymasspaper . I work with the familiar octet of light pseudo-scalar meson fields $`\pi _\alpha `$ with $`\alpha =1\mathrm{}8`$. In a standard way (see for example km ) I consider an effective field theory defined in terms of the SU(3) valued group element

$$\mathrm{\Sigma }=\mathrm{exp}(i\pi _\alpha \lambda _\alpha /f_\pi )SU(3).$$
(1)
Here the $`\lambda _\alpha `$ are the usual Gell-Mann matrices which generate the flavor group and $`f_\pi `$ is a dimensional constant with a phenomenological value of about 93 MeV. I follow the normalization convention that $`\mathrm{Tr}\lambda _\alpha \lambda _\beta =2\delta _{\alpha \beta }`$. The neutral pion and the eta meson play a special role in the later discussion; they are the coefficients of the commuting generators

$$\lambda _3=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)$$
(2)
and

$$\lambda _8=\frac{1}{\sqrt{3}}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 2\end{array}\right),$$
(3)
respectively. In the chiral limit of vanishing quark masses, we model the interactions of the eight massless Goldstone bosons with the effective Lagrangian density

$$L_0=\frac{f_\pi ^2}{4}\mathrm{Tr}(_\mu \mathrm{\Sigma }^{}_\mu \mathrm{\Sigma }).$$
(4)
The non-linear constraint of $`\mathrm{\Sigma }`$ onto the group SU(3) makes this theory non-renormalizable. It is to be understood only as the starting point for an expansion of particle interactions in powers of their masses and momenta. Expanding Eq. (4) to second order in the meson fields gives the conventional kinetic terms for our eight mesons.

This theory is invariant under parity and charge conjugation. These operators generate the transformations

$$\begin{array}{cc}& P:\mathrm{\Sigma }\mathrm{\Sigma }^1\\ & CP:\mathrm{\Sigma }\mathrm{\Sigma }^{}\end{array}$$
(5)
where the operation $``$ refers to complex conjugation. The eight meson fields are pseudo-scalars. The neutral pion and the eta meson are both even under charge conjugation.

With massless quarks, the underlying quark-gluon theory has a chiral symmetry under

$$\begin{array}{cc}& \psi _L\psi _Lg_L\\ & \psi _R\psi _Rg_R.\end{array}$$
(6)
Here $`(g_L,g_R)`$ is in $`(SU(3)\times SU(3))`$ and $`\psi _{L,R}`$ represent the chiral components of the quark fields, with flavor indices understood. This symmetry is expected to be broken spontaneously to a vector SU(3) via a vacuum expectation value for $`\overline{\psi }_L\psi _R`$. This motivates the sigma model through the identification

$$0|\overline{\psi }_L\psi _R|0v\mathrm{\Sigma }.$$
(7)
The quantity $`v`$, of dimension mass cubed, characterizes the strength of the spontaneous breaking of this symmetry. Thus our effective field transforms under the chiral symmetry as

$$\mathrm{\Sigma }g_L^{}\mathrm{\Sigma }g_R.$$
(8)
Our initial Lagrangian density is the simplest non-trivial expression invariant under this symmetry.

The quark masses break the chiral symmetry explicitly. From the analogy in Eq. (7), these are introduced through a 3 by 3 mass matrix $`M`$ appearing in a potential term added to the Lagrangian density

$$L=L_0v\mathrm{Re}\mathrm{Tr}(\mathrm{\Sigma }M).$$
(9)
Here $`v`$ is the same dimensionful factor appearing in Eq. (7). The chiral symmetry of our starting theory shows the physical equivalence of a given mass matrix $`M`$ with a rotated matrix $`g_R^{}Mg_L`$. Using this freedom we can put the mass matrix into a standard form. I will assume it is diagonal with increasing eigenvalues

$$M=\left(\begin{array}{ccc}m_u& 0& 0\\ 0& m_d& 0\\ 0& 0& m_s\end{array}\right)$$
(10)
representing the up, down, and strange quark masses. Note that this matrix has both singlet and octet parts under the vector flavor symmetry

$$M=\frac{m_u+m_d+m_s}{3}+\frac{m_um_d}{2}\lambda _3+\frac{m_u+m_d2m_s}{2\sqrt{3}}\lambda _8.$$
(11)
In general the mass matrix can still be complex. The chiral symmetry allows us to move phases between the masses, but the determinant of $`M`$ is invariant. Under charge conjugation the mass term would only be invariant if $`M=M^{}`$. If $`|M|`$ is not real, then its phase is the famous CP violating parameter usually associated with topological structure in the gauge fields. For the moment I take all quark masses as real. Since I am looking for spontaneous CP violation, I consider the case where there is no explicit CP violation.

To lowest order the masses of the pseudo-scalar mesons appear on expanding the mass term quadratically in the meson fields. This generates an effective mass matrix for the eight mesons

$$_{\alpha \beta }\mathrm{Re}\mathrm{Tr}\lambda _\alpha \lambda _\beta M.$$
(12)
The isospin-breaking up-down mass difference plays a crucial role in the later discussion. This gives this matrix an off diagonal piece mixing the $`\pi _0`$ and the $`\eta `$

$$_{3,8}m_um_d.$$
(13)
The eigenvalues of this matrix give the standard mass relations

$$\begin{array}{c}m_{\pi _0}^2\frac{2}{3}\left(m_u+m_d+m_s\sqrt{m_u^2+m_d^2+m_s^2m_um_dm_um_sm_dm_s}\right)\\ m_{\pi _+}^2=m_\pi _{}^2m_u+m_d\\ m_{K_+}^2=m_K_{}^2m_u+m_s\\ m_{K_0}^2=m_{\overline{K}_0}^2m_d+m_s\\ m_\eta ^2\frac{2}{3}\left(m_u+m_d+m_s+\sqrt{m_u^2+m_d^2+m_s^2m_um_dm_um_sm_dm_s}\right).\end{array}$$
(14)
Here I label the mesons with their conventional names.

Redundancies in these relations test the validity of the model. For example, comparing two expressions for the sum of the three quark masses

$$\frac{2(m_{\pi _+}^2+m_{K_+}^2+m_{K_0}^2)}{3(m_\eta ^2+m_{\pi _0}^2)}1.07$$
(15)
suggests the symmetry should be good to a few percent. Further ratios of meson masses then give estimates for the ratios of the quark masses km ; weinberg ; leutwyler . For one such combination, look at

$$\frac{m_u}{m_d}=\frac{m_{\pi ^+}^2+m_{K_+}^2m_{K_0}^2}{m_{\pi ^+}^2m_{K_+}^2+m_{K_0}^2}0.66.$$
(16)
This particular combination is polluted by electromagnetic effects; another combination partially cancels such while ignoring small $`m_um_d/m_s`$ corrections

$$\frac{m_u}{m_d}=\frac{2m_{\pi ^0}^2m_{\pi ^+}^2+m_{K_+}^2m_{K_0}^2}{m_{\pi ^+}^2m_{K_+}^2+m_{K_0}^2}0.55.$$
(17)
Later I will comment on a third combination for this ratio. For the strange quark, one can take

$$\frac{2m_s}{m_u+m_d}=\frac{m_{K_+}^2+m_{K_0}^2m_{\pi ^+}^2}{m_{\pi ^+}^2}26.$$
(18)
## III Spontaneous CP violation

So far all this is standard. Now I vary the quark masses and look for interesting phenomena. In particular, I want to find spontaneous breaking of the CP symmetry. Normally the $`\mathrm{\Sigma }`$ field fluctuates around the identity in SU(3). However, for some values of the quark masses this ceases to be true. When the vacuum expectation of $`\mathrm{\Sigma }`$ deviates from the identity, some of the meson fields acquire expectation values. As they are pseudo-scalars, this necessarily involves a breakdown of parity, as noted by Dashen dashen .

To explore this possibility, I concentrate on the lightest meson from Eq. (14), the $`\pi _0`$. From Eq. (14) we can calculate the product of the $`\pi _0`$ and $`\eta `$ masses

$$m_{\pi _0}^2m_\eta ^2m_um_d+m_um_s+m_dm_s.$$
(19)
Whenever

$$m_u=\frac{m_sm_d}{m_s+m_d}$$
(20)
the $`\pi _0`$ mass vanishes. For increasingly negative up-quark masses, our simple expansion around vanishing pseudo-scalar meson fields fails. The vacuum is no longer approximated by fluctuations of $`\mathrm{\Sigma }`$ around the unit matrix; instead it fluctuates about an SU(3) matrix of form

$$\mathrm{\Sigma }=\left(\begin{array}{ccc}e^{i\varphi _1}& 0& 0\\ 0& e^{i\varphi _2}& 0\\ 0& 0& e^{i\varphi _1i\varphi _2}\end{array}\right)$$
(21)
where the phases satisfy

$$m_u\mathrm{sin}(\varphi _1)=m_d\mathrm{sin}(\varphi _2)=m_s\mathrm{sin}(\varphi _1+\varphi _2).$$
(22)
There are two minimum action solutions, differing by flipping the signs of these angles. The transition is a continuous one, with $`\mathrm{\Sigma }`$ going smoothly to the identity as the boundary given by Eq. (20) is approached.

In the new vacuum the neutral pseudo-scalar meson fields acquire expectation values. As the neutral pion is CP odd, we spontaneously break this symmetry. This will have various experimental consequences, for example eta decay into two pions becomes allowed since a virtual third pion can be absorbed by the vacuum. Fig. (1) sketches the inferred phase diagram as a function of the up and down quark masses. Chiral rotations insure a symmetry under the flipping of the signs of both quark masses.

At first sight the appearance of the CP violating phase at negative up quark mass may seem surprising. Naively in perturbation theory the sign of a fermion mass can be rotated away by a redefinition $`\psi \gamma _5\psi `$. However this rotation is anomalous, making the sign of the quark mass observable. A more general complex phase in the mass would also have physical consequences, i.e. explicit CP violation. With real quark masses the underlying Lagrangian is CP invariant, but the above discussion shows that there exists a large region where the ground state spontaneously breaks this symmetry.

Vafa and Witten vafawitten  argued on rather general conditions that CP could not be spontaneously broken in the strong interactions. However their argument makes positivity assumptions on the path integral measure. When a quark mass is negative, the fermion determinant need not be positive for all gauge configurations; in this case their assumptions fail.

The possible existence of this phase was anticipated on the lattice some time ago by Aoki aoki . For the one flavor case he found this parity breaking phase with Wilson lattice gauge fermions. He went on to discuss also two flavors, finding both flavor and parity symmetry breaking. The latter case is now regarded as a lattice artifact of Wilson fermions. For a review of these issues see myreview .

In conventional discussions of CP non-invariance in the strong interactions witten  appears a phase $`e^{i\theta }`$ appearing on tunneling between topologically distinct gauge field configurations. The famous U(1) anomaly formally allows us to move this phase into the determinant of the quark mass matrix. After rotating all phases into the up-quark mass, we see that our spontaneous breaking of CP is occurring at an angle $`\theta =\pi `$.

A crucial observation is that when the down quark mass is positive, the CP violating phase does not appear for up-quark masses greater than a non-zero minimum value. There exists a finite gap with $`\theta =\pi `$ without this symmetry breaking. The chiral model predicts a smooth behavior as the up-quark mass passes through zero. This is the main motivation for the discussion in section V on the meaning of a vanishing up-quark mass. Indeed, from the effective Lagrangian point of view, the real and imaginary parts of the quark mass are independent parameters. The absence of experimental evidence for strong CP violation suggests that the imaginary part of the quark mass matrix vanishes, but says nothing about the real part.

An interesting special case occurs when the up and down quarks have the same magnitude but opposite sign for their masses, i.e. $`m_u=m_d`$. In this situation it is illuminating to rotate the minus sign into the phase of the strange quark. Then the up and down quark are degenerate, and we have restored an exact vector $`SU(2)`$ flavor symmetry. The excitation spectrum will show three degenerate pions, but they will not be massless due to what might be thought of a vacuum condensate of eta particles.

## IV Including the $`\eta ^{}`$

The above discussion was entirely in terms of the pseudo-scalar mesons that become Goldstone bosons in the chiral limit. One might wonder how higher states can influence this phase structure. Of particular concern is the $`\eta ^{}`$ meson associated with the anomalous $`U(1)`$ symmetry present in the classical quark-gluon Lagrangian. Non-perturbative processes, including topologically non-trivial gauge field configurations, are well known to generate a mass for this particle. I will now argue that, while this state can shift masses due to mixing with the lighter mesons, it does not make a qualitative difference in the existence of a phase with spontaneous CP violation.

The easiest way to introduce the $`\eta ^{}`$ into the effective theory is to promote the group element $`\mathrm{\Sigma }`$ to an element of $`U(3)`$ via an overall phase factor. Thus I generalize Eq. (1) to

$$\mathrm{\Sigma }=\mathrm{exp}\left(i\pi _\alpha \lambda _\alpha /f_\pi +i\sqrt{\frac{2}{3}}\eta ^{}/f_\pi \right)U(3).$$
(23)
The factor $`\sqrt{2/3}`$ gives the $`\eta ^{}`$ field the same normalization as the $`\pi `$ fields. Our starting kinetic Lagrangian in Eq. (4) would have this particle also be massless. One way to fix this deficiency is to mimic the anomaly with a term proportional to the determinant of $`\mathrm{\Sigma }`$

$$L_0=\frac{f_\pi ^2}{4}\mathrm{Tr}(_\mu \mathrm{\Sigma }^{}_\mu \mathrm{\Sigma })C|\mathrm{\Sigma }|.$$
(24)
The parameter $`C`$ parameterizes the strength of the anomaly in the $`U(1)`$ factor.

Now if we include the mass term exactly as before, additional mixing occurs between the $`\eta ^{}`$, the $`\pi _0`$, and the $`\eta `$. The corresponding mixing matrix takes the form

$$\left(\begin{array}{ccc}m_u+m_d& \frac{m_um_d}{\sqrt{3}}& \sqrt{\frac{2}{3}}(m_um_d)\\ \frac{m_um_d}{\sqrt{3}}& \frac{m_u+m_d+4m_s}{3}& \frac{\sqrt{2}(m_u+m_d2m_s)}{3}\\ \sqrt{\frac{2}{3}}(m_um_d)& \frac{\sqrt{2}(m_u+m_d2m_s)}{3}& \frac{2}{3}m_a\end{array}\right)$$
(25)
where $`m_a`$ characterizes the contribution of the non-perturbative physics to the $`\eta ^{}`$ mass. This should have a value of order the strong interaction scale; in particular, it should be large compared to at least the up and down quark masses. The two by two matrix in the upper left of this expression is exactly what is diagonalized to find the neutral pion and eta masses in Eq. (14).

The boundary of the CP violating phase occurs where the determinant of this matrix vanishes. This modifies Eq. (19) to

$$m_{\pi _0}^2m_\eta ^2m_\eta ^{}^2m_a(m_um_d+m_um_s+m_dm_s)m_u(m_dm_s)^2m_d(m_um_s)^2m_s(m_um_d)^2.$$
(26)
The boundary shifts slightly from the earlier result, but still passes through the origin, leaving Fig. (1) qualitatively unchanged.

## V Can the up quark be massless?

A oft proposed solution to the strong CP problem banks ; leutwyler0 ; cohen ; fleming  asks whether $`m_u=0`$. From the effective Lagrangian point of view, this appears to be an artificial setting of two parameters to zero, the real and imaginary parts of the quark mass. The earlier discussion shows that nothing special is expected to happen as the real part of the quark mass goes through zero. It is only the imaginary part that should vanish for CP to be a good symmetry, at least when the up-quark mass is larger than the value giving spontaneous breaking. But speculations on a vanishing up-quark mass continue, so it is interesting to ask if this can be given physical meaning. In this section I investigate precisely what is meant by a quark mass, and what $`m_u=0`$ would mean. I will conclude that the question of whether $`m_u`$ could vanish is ill posed. This is not relevant for most phenomenological purposes, but is unacceptable for solving something fundamental, like the strong CP problem. Ref. banks  raises some of these issues, pursuing $`m_u=0`$ anyway as an accidental symmetry.

The conventional description of this phenomenon changes variables from the complex quark mass to polar coordinates involving the magnitude of the up-quark mass and its phase. I argue below that non-perturbative effects give a scheme-dependent additive shift in the real part of the up-quark mass. Such a shift exposes the singular nature of this selection of coordinates. To emphasize the point, imagine changing variables to the magnitude and phase of $`(m_u4\mathrm{MeV})`$. If the up-quark mass is 4 MeV, then in such coordinates the theory no longer depends on the phase. From this point of view $`m_u=4\mathrm{MeV}`$ is an equally good solution to the strong CP issue as is $`m_u=0`$.

While phenomenology, i.e. Eq. (17), seems to suggest that the up quark is not massless, there remains considerable freedom in extracting that ratio from the masses of the pseudo-scalar mesons. From Eq. (14), the sum of the squares of the $`\eta `$ and $`\pi _0`$ masses should be proportional to the sum of the three quark masses. Subtracting off the neutral kaon mass should leave just the up quark. Thus motivated, look at

$$\frac{m_u}{m_d}=\frac{3(m_\eta ^2+m_{\pi _0}^2)/22m_{K_0}^2}{m_{\pi ^+}^2m_{K_+}^2+m_{K_0}^2}0.8.$$
(27)
Thus even the sign of the up-quark mass is ambiguous. This example is perhaps a bit extreme since it ignores the influence of mixing with the eta prime in the numerator. But it shows that there does exist a substantial phenomenological uncertainty in the quark masses. More formal attempts to extend the naive quark mass ratio estimates to higher orders in the chiral expansion have shown fundamental ambiguities in the definition of the quark masses km .

If two quark masses were to vanish simultaneously, then we would have exactly massless pions, Goldstone bosons for the resulting flavored chiral symmetry. In this case the concept of zero quark mass has definite physical consequences. But here I concentrate on whether the concept of only a single massless quark has any meaning. While I could carry along the baggage of the heavier quarks, let me simplify the discussion and consider the theory reduced to a single flavor of quark.

As stated in the introduction, because renormalization is required, the continuum concept of an “underlying basic Lagrangian” does not exist. The continuum theory is specified in terms of basic symmetries and a few renormalized parameters. Because of anomalies, a single massless quark does not correspond to any symmetry. In practice the definition of a field theory relies on a limiting process from a cutoff version. As the lattice is the only well understood non-perturbative cutoff, it provides the most natural framework for such a definition. But any regulator must accommodate the known chiral anomalies, and thus there must be chiral symmetry breaking terms in the cutoff theory. These chiral breaking effects come in many guises. With a Pauli-Villars scheme, there is a heavy regulator field. With dimensional regularization the anomaly is hidden in the fermionic measure. For Wilson lattice gauge theory there is the famous Wilson term. With domain wall fermions there is a residual mass from a finite fifth dimension. With overlap fermions things are hidden in a combination of the measure and a certain non-uniqueness of the operator. I will return to this last case shortly.

The lattice regulator involves introducing a dimensionful parameter, the lattice spacing $`a`$. This feature is not special to the lattice. The scale anomaly is responsible for masses of hadrons such as the proton and glueballs, even in the massless quark limit. For such physics, any complete regulator must introduce a scale.

The renormalization process tunes all relevant bare parameters as a function of the cutoff to fix a set of renormalized quantities. In the case of the strong interactions, the bare gauge coupling is driven to zero by asymptotic freedom. Its cutoff dependence is absorbed into an overall scale in a well known way via the phenomenon of dimensional transmutation cw . The only other parameters of the strong interactions are the quark masses. For these one inputs a few particle masses to finally determine the continuum theory uniquely. For the three flavor theory the most natural observables to fix these parameters are the masses of the pseudo-scalar mesons.

In the one flavor theory there are no Goldstone bosons, but massive mesons and baryons should exist. I need some physical parameter with which to carry out the renormalization of the quark mass. For this purpose I choose the ratio of the lightest boson mass to the lightest baryon mass. As both are expected to be stable, this precludes any ambiguity from particle widths. Calling the lightest boson the $`\eta `$ and the baryon $`p`$, I define

$$r=\frac{m_\eta }{m_p}.$$
(28)
I expect to be able to adjust this parameter via the quark mass, which should be tuned to give the desired value. It should be possible to give this ratio any value throughout the range from $`r=0`$ at the boundary of the above CP violating phase to $`2/3`$ in the heavy quark limit.

With a cutoff in place, I can in principle determine this ratio given any values for the bare quark mass and bare coupling. For pedagogy, let me assume a lattice regulator and trade these parameters for the lattice spacing $`a`$ and the quark mass in lattice units, $`m_qa`$. Both of these quantities go to zero in the continuum limit. The renormalization prescription is to select a desired value of $`r`$ and follow the contour with this value in the $`(a,m_q)`$ plane towards the origin. This process is sketched in Fig. (2). Perturbative divergences in the bare quark mass appear in the fact that these contours approach the origin with zero slope.

When the cutoff is in place, we expect prescription dependent artifacts. In particular, the precise locations of the constant $`r`$ contours will depend on details of the formulation. Holding the bare quark mass at zero will cross a variety of $`r`$ contours, with none obviously favored as the origin is approached. Different cutoff schemes will give different continuum limits for $`m_u=0`$. Alternatively, $`m_u=0`$ in one scheme could give the same continuum limit as some other scheme with a non-vanishing $`m_u`$. A non-perturbative additive shift in the up-quark mass distinguishes the two schemes. Asking that the up-quark mass vanishes is unphysical.

With two or more degenerate flavors there will be one special contour where the lightest meson does represent a Goldstone boson. With the Wilson fermion formulation, the quark mass axis is represented by the hopping parameter. As this particular cutoff explicitly breaks chiral symmetry, the critical hopping parameter, where the meson mass vanishes, is renormalized away from its value in the continuum limit.

Recently there has been considerable progress with lattice fermion formulations that preserve a remnant of exact chiral symmetry myreview . With such, the two flavor theory will have the $`r=0`$ contour naturally preserved as the $`m_q=0`$ axis. The crucial point is that this is not true for the one flavor theory. Dynamics generates a mass for the pseudo-scalar meson; thus, the $`m_q=0`$ axis will cut through various finite values of $`r`$. An interesting question is whether, as we take the lattice spacing to zero along this axis, some physical value of $`r`$ will be picked out as special and corresponding to vanishing quark mass. That this is unlikely follows from the non-uniqueness of these chiral lattice operators. For example, the overlap operator overlap  is constructed by a projection process from the conventional Wilson lattice operator. The latter has a mass parameter which is to be chosen in a particular domain. On changing this parameter, the massless Dirac operator still satisfies the Ginsparg-Wilson relation gw , but this condition does not guarantee that the contours of constant $`r`$ in Fig. (2) will not shift. Thus the horizontal axis is not expected to select one contour as special. Again, holding $`m_u=0`$ is not expected to give a unique continuum theory.

To see that this non-universality is indeed expected, consider that the dynamics of the one flavor case generates a mass gap in the $`\eta `$ channel. This means that the eigenvalues of the Dirac operator important to low energy physics are are not near the origin, but dynamically driven a finite distance away. Changing the projection procedure to generate the overlap operator will modify the size of this gap, changing the $`\eta `$ mass and the ratio $`r`$.

One might attempt to define $`m_u=0`$ as the point where the topological susceptibility vanishes. For this purpose, the overlap operator provides a natural definition of the gauge field topology. With this prescription a vanishing topological susceptibility is synonymous with vanishing quark mass for that particular overlap operator. As the result depends on the operator chosen, with one flavor the susceptibility is no more physical than the quark mass. As a rather abstract concept, the topological susceptibility is not directly measurable in scattering processes for physical particles.

This non-perturbative ambiguity in the quark mass carries over to the explicitly CP violating case where the mass is complex with phase $`\theta `$. The above discussion shows that even the sign of the up quark mass is ambiguous, indicating that different schemes can have an ambiguity of $`\theta `$ between 0 and $`\pi `$, a particularly severe example. For other values of $`\theta `$, to fix the continuum theory uniquely we need to introduce another renormalized quantity. For example, this could be a three meson coupling or the electric dipole moment of a baryon. Beyond the lowest order in the chiral expansion, the precise dependence of the renormalized parameter on $`\theta `$ is scheme-dependent.

While the use of the lattice provides a framework for precise discussion, the source of the additive shift in the up-quark mass follows qualitatively from non-perturbative classical gauge configurations, i.e. “pseudoparticles” or “instantons” mcarthurgeorgi . As shown some time ago by ’t Hooft thooft , these configurations generate an effective multi-fermion vertex where all flavors of quark flip their spin. If we take this vertex and tie together the massive quark lines with mass terms, then the resulting process generates an effective mass term for the light quark. The strength of this term is proportional to the product of the masses of the more massive quarks. This process is illustrated in Fig. (3). As the strength of the effective interaction is scheme and scale dependent, so is the resulting light quark mass.

## VI Final remarks

While I have been exploring rather unphysical regions in parameter space, these observations do raise some issues for practical lattice calculations of hadronic physics. Current simulations are done at relatively heavy values for the quark masses. This is because the known fermion algorithms tend to converge rather slowly at light quark masses. Extrapolations by several tens of MeV are needed to reach physical quark masses, and these extrapolations tend to be made in the context of chiral perturbation theory. The presence of a CP violating phase quite near the physical values for the quark masses indicates a strong variation in the vacuum state with a rather small change in the up-quark mass; indeed, less than a 10 MeV change in the traditionally determined up-quark mass can drastically change the low energy spectrum. Most simulations consider degenerate quarks, and chiral extrapolations so far have been quite successful. But some quantities, namely certain baryonic properties thomas , do seem to require rather strong variations as the chiral limit is approached. These effects and the strong dependence on the up-quark mass may be related.

Another issue is the validity of current simulation algorithms with non-degenerate quarks. With an even number of degenerate flavors the fermion determinant is positive and can contribute to a measure for Monte Carlo simulations. With light non-degenerate quarks the positivity of this determinant is not guaranteed. Indeed, the CP violation can occur only when the fermions contribute large phases to the path integral. Current algorithms for dealing with non-degenerate quarks milc  take a root of the determinant with multiple flavors. In this process any possible phases are ignored. Such an algorithm is incapable of seeing any of the CP violating phenomena discussed here. This point may not be too serious in practice since the up and down quarks are nearly degenerate and the strange quark is fairly heavy. But these issues should serve as a warning that things might not work as well as we want.

## Acknowledgements

I have benefited from lively discussions with many colleagues, in particular F. Berruto, T. Blum, H. Neuberger and G. Senjanovic. This manuscript has been authored under contract number DE-AC02-98CH10886 with the U.S. Department of Energy. Accordingly, the U.S. Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes.