text:
# Channel Formation and Intermediate Range Order in Sodium Silicate Melts and Glasses

## Abstract

We use inelastic neutron scattering and molecular dynamics simulation to investigate the interplay between the structure and the fast sodium ion diffusion in various sodium silicates. With increasing temperature and decreasing density the structure factors exhibit an emerging prepeak around 0.9 Å-1. We show, that this prepeak has its origin in the formation of sodium rich channels in the static structure. The channels serve as preferential ion conducting pathways in the relative immobile Si–O matrix. On cooling below the glass transition this intermediate range order is frozen in.

Alkali silicates are a paradigm for multicomponent glass forming systems in which even at typical melt temperatures the mobility of the alkali ions exceeds that of the Si–O network by orders of magnitude JoBB51 ; GuKi66 ; BrFr88 . Although such a high mobility of the alkali ions should be reflected in the underlying structure, up to now, it is still unclear how the alkali atoms are built into the Si–O network. A distribution where an increasing alkali content causes increasing homogenous disruptions of the SiO4 tetrahedal network of pure silica Zac32 ; WaBi38  does not agree with the observed drastic difference in mobility. In addition, this distribution is in conflict with the highly non–linear dependence of the viscosity on alkali concentration MaSS83 ; KnDS94 . These observations lead to ideas that propose the existence of preferential ion conducting pathways in a Si–O matrix AnCT82 ; Gre85 .

Here, we present evidence for the existence of sodium diffusion channels in the static structure. Using neutron scattering experiments and molecular dynamics simulations on various sodium silicate melts and glasses, we investigate the structure that provides the mobility of the alkali ions. We demonstrate the accuracy of the molecular dynamics model by our neutron data. Then we make use of the detailed information obtained by molecular dynamics to elucidate the observed intermediate range order as seen by neutron scattering.

Sodium di- (NS2), tri- (NS3) and tetra- (NS4) silicate glasses were synthesized from ultrapure Na2CO3 and SiO2 powders by fusion at 1500 K. For the inelastic neutron scattering experiments glassy samples were encapsulated in Pt cells giving an annular sample geometry of 40 mm in height, 22.5 mm in diameter and a 1.25 mm wall thickness. Measurements were performed at room temperature and on liquid samples at temperatures between 1200 K and 1600 K on the neutron time–of–flight spectrometer IN 6 of the Institut Laue-Langevin in Grenoble. An incident neutron wavelength of $`\lambda =5.9`$ Å-1 yielded an energy resolution of $`\delta E=50\mu \text{eV}`$ (FWHM) and an accessible momentum transfer at zero energy transfer of $`q=0.21.75`$ Å-1. The raw data reduction consisted of a normalization to a vanadium standard and a correction for self absorption and container scattering.

At 1600 K the structural relaxation of the Si–O network for the investigated sodium silicates is on a nanosecond scale MeSD02 ; KnDS94  – too slow to be resolved on IN 6. Therefore, the dynamic structure factor $`S(q,\omega )`$, depending on momentum $`q`$ and energy transfer $`\mathrm{}\omega `$, exhibits a strong elastic line with an elastic structure factor $`S(q,\omega =0)`$ that is in good approximation the static structure factor times the Debye–Waller factor $`f(q)`$. Whereas scattering on the Si and O atoms is exclusively coherent, scattering on sodium is coherent and incoherent. The incoherent scattering reflects itself in a flat background in $`S(q,\omega =0)`$.

The potential that was used for the molecular dynamics simulations is a slight modification of the pair potential by Kramer et al. KrMS91  which is based on ab initio calculations. More details on this potential can be found in Ref. HoKB99 . The simulations were done at various densities for systems of $`N=8064`$ particles ($`N_{\mathrm{Si}}=2016`$, $`N_{\mathrm{Na}}=1344`$, and $`N_\mathrm{O}=4704`$). At the densities $`2.2`$ g/cm3 and $`2.37`$ g/cm3 systems at the temperature $`T=2100`$ K, i.e. in the liquid state, were equilibrated for 3.3 ns followed by microcanonical production runs of the same length (using the velocity form of the Verlet algorithm with a time step of 1.6 fs). In each case two independent runs were done in order to improve the statistics. The equilibrated systems at the density $`\rho =2.37`$ g/cm3 were also quenched to a glass state at $`T=300`$ K using a cooling rate of about $`10^{12}`$ K/s. In the latter simulation the pressure was kept constant at ambient pressure which yielded a density of $`2.46`$ g/cm3 at 300 K. Note that the densities $`\rho =2.2`$ g/cm3 and $`\rho =2.46`$ g/cm3 are close to the experimental densities at $`T=2100`$ K and $`T=300`$ K, respectively MaSS83 .

The central quantity for the determination of the elastic structure factor from the simulation is the time–dependent dynamic structure factor defined by

$$S(q,t)=\frac{N}{_\alpha N_\alpha b_\alpha ^2}\underset{\alpha ,\beta }{}b_\alpha b_\beta S_{\alpha \beta }(q,t)\alpha ,\beta \{\mathrm{Si},\mathrm{Na},\mathrm{O}\}$$
(1)
with

$$S_{\alpha \beta }(q,t)=\frac{1}{N}\underset{k=1}{\overset{N_\alpha }{}}\underset{l=1}{\overset{N_\beta }{}}\mathrm{exp}[i𝐪(𝐫_k(t)𝐫_l(0))],$$
(2)
$`𝐫_i(t)`$ being the position vector of particle $`i`$ at time $`t`$. In Eq. (1) $`b_\alpha `$ denotes the experimental coherent scattering lengths of particle species $`\alpha `$. The elastic structure factors were estimated by the equation $`S(q,\omega =0)=f(q)S(q)`$ (with $`S(q)S(q,t=0)`$) where the Debye–Waller factor $`f(q)`$ was determined by fitting in the $`\beta `$relaxation regime the normalized quantity $`F(q,t)=S(q,t)/S(q)`$ to a von Schweidler law (for details see Ref. HoK02 ). Note that we have not taken into account the incoherent scattering on sodium in the calculation of $`S(q,\omega =0)`$ because this contribution is not expected to change the shape of the function $`S(q,\omega =0)`$ significantly.

Fig. 1a displays the elastic structure factor of glassy and viscous sodium trisilicate as seen by inelastic neutron scattering. The maximum at $`1.7`$Å-1 corresponds to the disrupted tetrahedal Si–O network structure of alkali silicate WaSu77 ; StMD95 . With increasing temperature the elastic scattered intensity at this wave-vector is decreasing at the expense of an increasing inelastic scattering. In contrast $`S(q,\omega =0)`$ displays an additional prepeak at $`0.9`$Å-1 which grows with increasing temperature, indicating an enhanced probability for a correlation of atomic arrangements at distances around 6-8 Å.

In sodium disilicate MeSD02  and sodium tetrasilicate neutron scattering reveals a similar behaviour: The elastic structure factors exhibit a growing prepeak at $`0.9`$ Å-1 on temperature increase. With increasing sodium content the height of the prepeak is increasing, whereas its position remains fairly unaffected by this large change in composition (Fig. 1b). We note, that for a monoatomic system a change in the structure that is due to thermal expansion would lead – beside a slight shift of the peak positions towards a smaller $`q`$ value – to a decreasing intensity of the peak height with increasing temperature. Thus, the emerging prepeak in the elastic structure factors as shown in Fig. 1 does not seem to be merely caused by thermal expansion but might be related to a change in the underlying structure. In the following, we clarify this behaviour of the elastic structure factors by means of our molecular dynamics simulations.

Fig. 2 displays the “simulated” neutron scattering elastic structure factor of sodium trisilicate. At the experimental densities, $`S(q,\omega =0)`$ exhibits the same behavior as in the neutron scattering experiment incsca : Whereas around 1.7 Å-1 the elastic scattered intensity decreases with increasing temperature (as expected), around 0.9 Å-1 a shoulder is present at $`T=2100`$ K which is nearly absent in the glass at $`T=300`$ K. Also shown in Fig. 2 is the elastic structure factor at the density of $`2.37`$ g/cm3 at $`T=2100`$ K. Here the shoulder around 0.9 Å-1 has a smaller amplitude than in the corresponding case at the experimental density. This gives evidence that the possible structural changes leading to a more pronounced appearance of the feature around 0.9 Å-1 in $`S(q,\omega =0)`$ are not due to a change in temperature but due to the change in the density. Thus, the question arises whether the shoulder around 0.9 Å-1 is related to structural features that tend to disappear with increasing density.

The partial static structure factors $`S_{\alpha \beta }(q)`$ at the experimental densities for different correlations are shown in Fig. 3: A well–pronounced peak is present in $`S_{\mathrm{Si},\mathrm{Na}}(q)`$ and $`S_{\mathrm{Na},\mathrm{Na}}(q)`$ at $`T=2100`$ K and at $`T=300`$ K. Thus, the structure that leads to the peak at $`q_1=0.9`$ Å-1 does not at all disappear with increasing density. That the feature at $`q_1`$ seems to be absent in $`S(q,\omega =0)`$ at 300 K, is due to the fact that $`S(q,\omega =0)`$ is a linear combination of six different partial structure factors: On the one hand positive and negative contributions cancel each other (note the negative amplitude of $`S_{\mathrm{Si},\mathrm{Na}}(q)`$ at $`q_1`$), and on the other hand $`S_{\mathrm{O},\mathrm{O}}(q)`$ exhibits only a shoulder at $`q_1`$ which is less pronounced at $`T=300`$ K. Because oxygen is the majority component (about 60% of the particles in NS3) and the coherent scattering length of oxygen is significantly larger than that of sodium and silicon ($`b_\mathrm{O}/b_{\mathrm{Si}}1.4`$ and $`b_\mathrm{O}/b_{\mathrm{Na}}1.6`$), $`S_{\mathrm{O},\mathrm{O}}(q)`$ gives the major contribution to $`S(q,\omega =0)`$ and thus, the changes in $`S_{\mathrm{O},\mathrm{O}}(q)`$ are the main cause for the emerging prepeak in the experimental and simulated $`S(q,\omega =0)`$ with increasing temperature and decreasing density, respectively.

The feature at $`q_1`$ corresponds to a distribution of sodium ions that is not homogenous on a length scale of 6–8 Å. This is illustrated by the snapshot, Fig. 4, which shows a molecular dynamics configuration of NS3 at $`T=2100`$ K at the density 2.2 g/cm3. Here, the sodium ions are represented by blue spheres that are connected to each other by surfaces. For clarity, the Si and O atoms are shown as small yellow and red spheres, respectively, that are connected to each other by covalent bonds shown as sticks. Note that the size of the latter spheres does not mirror the real size of Si and O. The snapshot gives a clear picture of a network of sodium channels that percolate through the Si–O structure. The peak at $`q_1`$ in the static structure factor marks the characteristic length scale of the latter network of channels.

Recently, Jund et al. JuKJ01  have shown in a molecular dynamics simulation of NS4 that the sodium trajectories form a network of pockets and channels where the distance between the pockets is of the order of 5–8 Å. Then, Horbach et al. HoKB02  have given evidence in the case of NS2 that the latter network is reflected by the prepeak at $`q_1`$ in the static structure factor. Fig. 4 shows that the “structure of the sodium trajectories” is not only seen in the dynamics as claimed by Jund et al. JuKJ01 . The sodium trajectories are given by paths in a quasi–static structure of sodium diffusion channels, the lifetime of which is given by the characteristic relaxation time of the Si–O network.

A change of the Na2O content (from NS4 to NS2) or of the density does not significantly affect the position of the prepeak at $`q_1`$ as seen by molecular dynamics simulations and neutron scattering. The effect that the prepeak at $`q_1`$ becomes more pronounced in $`S(q,\omega =0)`$ with increasing Na2O concentration (Fig. 1b) can be also explained by the simulation HoKB99 : The amplitude of the peak around $`1.7`$ Å-1 decreases with an increasing sodium content indicating a stronger disruption of the Si–O network if more sodium ions are added to the system. Moreover, Na–Na correlations contribute with a larger weight in the total structure factor and thus also the amplitude of the prepeak at $`q_1`$ increases. These findings from the simulation explain the behavior of $`S(q,\omega =0)`$ as seen in Fig. 1b.

However, a change in the composition or the density is accompanied by a significant change in the disrupted Si–O network structure. This is indicated by the behaviour of the Si–O–Si angle: In the NS3 simulation at the experimental densities this angle increases from 134.5 at $`T=300`$ K (i.e. at the density of 2.46 g/cm3) to 143 (i.e. at the density of 2.2 g/cm3). Thus, the sodium diffusion channels are sustained over a broad range in temperature and sodium concentration because the packing of the disrupted Si–O tetrahedra is rather “flexible” with respect to the Si–O–Si angle. Of course, there should be a lower bound of sodium concentrations below which the latter rearrangements in the Si–O network are no longer possible. Indeed, at very low sodium oxide content (below about 5 mol%), a miscibility gap of Na2O and SiO2 even at typical melt temperatures has been reported mazurin .

Very recently, Lammert et al. LKH03  have investigated a lithium silicate melt by molecular dynamics simulation. They have shown that the diffusion mechanism for lithium is a kind of vacancy diffusion. Movies that we have done from the NS3 configurations such as Fig. 4 reveal a similar sodium diffusion mechanism. Thus, our picture of diffusion channels is consistent with the vacancy diffusion on a local length scale as proposed by Lammert et al. provided that also in lithium silicates diffusion channels exist. Indeed, our recent results of inelastic neutron scattering measurements KaMe03  and molecular dynamics simulations KnHB03  on various lithium and potassium bearing silicates indicate that preferential ion conduction pathways are a more general feature for alkali silicates. A neutron scattering experiment on a calcium silicate glass using Ca isotope substitution revealed a prepeak in the Ca–Ca partial structure factor GaEB91  at $`1.3`$ Å-1. In contrast to our findings in alkali silicates, the signal has been interpreted in terms of a glass structure based on ordered, densily packed domains. It remains to be seen, to what extent the existence of preferential ion conducting pathways in the static structure shown here for sodium silicates applies to other glass forming, oxide ion conductors.

In conclusion, inelastic neutron scattering experiments and molecular dynamics simulations on sodium silicate melts and glasses revealed the existence of sodium diffusion channels in the static structure. We find, that the formation of preferential ion conducting pathways is fairly unaffected by the sodium oxide concentration (in the studied composition range from sodium tetrasilicate to sodium disilicate), by a temperature change corresponding to the glass at 300 K up to the equilibrium melt, or a change of the density in the melt of about 8 %. In sodium silicates the values of the sodium ion diffusion at typical melt temperatures JoBB51 ; GuKi66 ; BrFr88  and the values of the viscosities MaSS83 ; KnDS94  depend relatively weakly on the sodium oxide concentration. The formation of sodium rich channels gives an explanation on a microscopic level for the observed macroscopic properties of mass transport in sodium silicates.

We thank Kurt Binder, Donald B. Dingwell, and Winfried Petry for their support and fruitful discussions. We are grateful to Hans Knoth for preparing Fig. 4. One of us (J. H.) was supported by the Deutsche Forschungsgemeinschaft (DFG) grant No HO2231/2–1.


format: text

text:
# Diffusion-based method for producing density equalizing maps

## 1 Introduction

Suppose we wish to represent on a map some data concerning, to take the most common example, the human population. For instance, we might wish to show votes in an election, incidence of a disease, number of cars, televisions, or phones in use, numbers of people falling in one group or another of the population, by age or income, or any of very many other variables of statistical, medical, or demographic interest. The typical course under such circumstances would be to choose one of the standard projections for the area of interest and plot the data on it using some color code or similar representation. Such maps however can be highly misleading. A plot of disease incidence, for example, will inevitably show high incidence in cities and low incidence in rural areas, solely because there are more people living in cities.

The obvious cure for this problem is to plot a fractional measure rather than raw incidence data; we plot some measure of the number of cases per capita, binned in segments small enough to give good spatial resolution but large enough to give reliable sampling. This however has its own problems, since it discards all information about where most of the cases are occurring. One case in a thousand means something entirely different in Sydney from what it means in Siberia.

What we would like is some representation of the data that factors out variations in the population density but, at the same time, shows how many cases are occurring in each region. It appears at first that these two goals are irreconcilable, but this is not the case. On a normal area-preserving or approximately area-preserving projection, such as a Mercator or Robinson projection, they are indeed irreconcilable. However, if we can construct a projection in which areas on the map are proportional not to areas on the ground but instead to human population (or whatever the corresponding variable might be for the problem in hand), then we can have our cake and eat it. Disease cases or other similar data plotted on such a projection will have the same density in areas with equal per capita incidence regardless of the population, since both the raw incidence rate and the area will scale with the population. However, each case or group of cases can still be represented individually, so it will be clear to the eye where most of the cases occur. Projections of this kind are known as value-by-area maps, density equalizing maps, or cartograms.

The construction of cartograms turns out to be a challenging undertaking. A variety of methods have been put forward, but none is entirely satisfactory. In particular, they often produce highly distorted maps that are difficult to read or projections that are badly behaved under some circumstances, with overlapping regions or strong dependence on coordinate axes. In many cases the methods proposed are also computationally demanding, sometimes taking hours to produce a single map. In this paper we propose a new method which is, we believe, intuitive, while also producing elegant, well-behaved, and useful cartograms, whose calculation makes relatively low demands on our computational resources.

## 2 Previous methods for constructing cartograms

Before introducing our own method for constructing cartograms, we outline in this section some of the previous methods that have been proposed. We will in most of the paper be considering cartograms based on population density, although the developments are easily generalized to other quantities. (We give one different example in Sec. 5.) The construction of a (flat two-dimensional) cartogram involves finding a transformation $`𝐫𝐓(𝐫)`$ of a plane to another plane such that the Jacobian $`(T_x,T_y)/(x,y)`$ of the transformation is proportional to some specified (population) density $`\rho (𝐫)`$ thus:

$$\frac{(T_x,T_y)}{(x,y)}\frac{T_x}{x}\frac{T_y}{y}\frac{T_x}{y}\frac{T_y}{x}=\frac{\rho (𝐫)}{\overline{\rho }},$$
(1)
where $`\overline{\rho }`$ is the mean population density averaged over the area $`A`$ to be mapped. (This choice of normalization for the Jacobian ensures that the total area before and after the transformation is the same.)

Eq. (1) does not determine the cartogram projection uniquely. To do that we need one more constraint: two constraints are needed to fix the projection for a two-dimensional cartogram. Different choices of the second constraint give different projections, and no single choice appears to be the obvious candidate, which is why many methods of making cartograms have been suggested. One idea is to demand conformal invariance under the cartogram transformation, i.e., to demand that angles be preserved locally. This is equivalent to demanding that the Cauchy–Riemann equations be satisfied, but this imposes two, not one, additional constraints and hence it is normally not possible to construct a conformally invariant cartogram.

In an attempt at least to minimize the distortion of angles, Tobler  took the first steps in the automated computer generation of cartograms in the late 1960s. He proposed a method in which the initial map is divided into small rectangular or hexagonal cells, each of which is then independently dilated or shrunk to a size proportional to its population content. Since each cell is scaled separately, the corners of adjacent cells don’t match afterwards. To re-establish a match, Tobler’s method takes a vector average over the positions of corresponding corners and draws a new map with the resulting distorted cells. The process is iterated until a fixed point of the transformation is reached. Although the principle is simple and intuitive it runs into practical problems. First, convergence tends to be rather slow because a node a few cells away from a population center will feel the effect of that center only after several iterations. Second, under some circumstances the transformation can produce overlapping or “folded” regions of the map, thereby ruining the topology. This problem can be corrected by introducing additional constraints but the result is a more complex algorithm with even slower run-times.

To increase the speed of the calculations, Dougenik et al.  introduced an algorithm where the borders of a cell move in response not only to local space requirements but also to “forces” exerted by other cells. Cells create force fields that diminish with distance from the cell and that are larger for cells that contain larger populations. These forces “push” other cells away from areas of high population in a manner reminiscent of the behavior of charged objects in electrostatics (although the authors do not use this metaphor). Again the positions are relaxed iteratively to achieve the final cartogram, and convergence is substantially faster than Tobler’s algorithm, although topological errors still cannot be ruled out.

Gusein-Zade and Tikunov  suggested a further twist that does away with the cells altogether and uses a continuous “displacement field” that measures the displacement of each point in the map. Areas of high population exert a repulsive force on this displacement field and the authors are able to derive a differential equation for the field, which they integrate numerically. The method is slow and somewhat arcane but produces some of the most attractive cartograms among the existing algorithms.

A more intuitive method, which makes use of cellular automata, has recently been proposed by Dorling . In this method the original map is drawn on a fine grid. On each iteration of the algorithm, cells lying on or close to the boundaries of regions are identified and if a neighboring region needs extra area those cells are reassigned to the neighbor. The procedure is iterated and the regions with greatest population grow slowly larger until an equilibrium is reached and no further changes are needed. The procedure is elegant and simple, but in practice it can distort shapes quite badly. One can add additional constraints on the shapes to make the maps more readable, but then the method quickly loses its main advantage, namely its simplicity.

Researchers have also experimented with several other methods. Kocmoud , for example, uses a mass-and-spring model acting on a map expressed as points and lines, with constraints applied to maintain certain topographic features such as angles or lengths. Due to its complexity, however, this algorithm is prohibitively slow. The method of Keim et al. , by contrast, is very fast, but achieves its speed primarily by working with polygonal maps that have been heavily simplified before beginning the computations, which unfortunately dispenses with many useful cartographic details. Finally, if one is willing to live with a non-contiguous cartogram—one in which regions adjacent in real life are not adjacent on the cartogram—then several quite simple methods give good results, such as Dorling’s circular cartograms . Other reviews and discussions of cartogram methods can be found in Refs. .

## 3 The diffusion cartogram

In this paper we propose a new type of cartogram, which might be called a “diffusion cartogram,” for reasons we now describe. It is a trivial observation that on a true population cartogram the population is necessarily uniform: once the areas of regions have been scaled to be proportional to their population then, by definition, population *density* is the same everywhere—hence the name “density-equalizing map.” Thus, one way to create a cartogram given a particular population density, is to allow population somehow to “flow away” from high-density areas into low-density ones, until the density is equalized everywhere. There is an obvious candidate process that achieves this, the linear diffusion process of elementary physics , and this is the basis of our method. We describe the population by a density function $`\rho (𝐫)`$, where $`𝐫`$ represents geographic position, and then allow this density to diffuse. In the limit of long time $`t\mathrm{}`$, the population density becomes uniform and so comes to rest, and its total displacement from start to finish determines the projection of the map necessary to produce a perfectly density-equalizing cartogram.

In the standard treatment of diffusion, the current density is given by

$$𝐉=𝐯(𝐫,t)\rho (𝐫,t),$$
(2)
where $`𝐯(𝐫,t)`$ and $`\rho (𝐫,t)`$ are the velocity and density respectively at position $`𝐫`$ and time $`t`$. Diffusion follows the gradient of the density field thus:

$$𝐉=\rho ,$$
(3)
meaning that the flow is always directed from regions of high density to regions of low and that the flow will be faster when the gradient is steeper. There is conventionally a diffusion constant in Eq. (3), but we can set this constant to 1 without affecting our results.

The diffusing population is also conserved locally so that

$$𝐉+\frac{\rho }{t}=0.$$
(4)
Combining Eqs. (2), (3), and (4) we then arrive at the familiar diffusion equation:

$$^2\rho \frac{\rho }{t}=0$$
(5)
and the expression for the velocity field in terms of the population density:

$$𝐯(𝐫,t)=\frac{\rho }{\rho }.$$
(6)
The calculation of the cartogram involves solving Eq. (5) for $`\rho (𝐫,t)`$ starting from the initial condition in which $`\rho `$ is equal to the given population density of the region of interest and then calculating the corresponding velocity field from Eq. (6). The cumulative displacement $`𝐫(t)`$ of any point on the map at time $`t`$ can be calculated by integrating the velocity field thus:

$$𝐫(t)=𝐫(0)+_0^t𝐯(𝐫,t^{})dt^{}.$$
(7)
In the limit $`t\mathrm{}`$ the set of such displacements for all points on the original map defines the cartogram.

Most of the time we are not interested in mapping the entire globe, but only some part of it, which means that the area of interest will have boundaries beyond which we don’t know or don’t care about the population density. These boundaries might be country borders or other political boundaries or coast lines. It would be inappropriate to represent regions outside the boundaries as having zero population, even if they are, like the ocean for instance, unpopulated, since this would cause arbitrary expansion of the cartogram as the population diffused into its uninhabited surroundings. (This is true of essentially all methods for constructing cartograms.) Instead, therefore, we apply a “neutral buoyancy” condition, floating the area of interest in a “sea” of uniform population density equal to the mean density of the map as a whole. This keeps the total area under consideration constant during the diffusion process.

The whole simulation, including the sea, is then enclosed in a box. For simplicity in this paper we will consider only rectangular boxes, as most others have also done. Provided the dimensions $`L_x`$ and $`L_y`$ of the box are substantially larger than the area to be mapped, the dimensions themselves do not matter. In the limit $`L_x,L_y\mathrm{}`$ the cartogram will be a unique deterministic mapping, independent of the coordinate system used, with no overlapping regions. In practice, we find that quite moderate system sizes are adequate—dimensions three to four times the linear extent of the area to be mapped appear to give good results.

We also need to choose boundary conditions on the walls of the box. These too have no great effect on the results, provided the size of the box is reasonably generous, and we have found a good choice to be the Neumann boundary conditions in which there is no flow of population through the walls of the box.

The considerations above completely specify our method and are intuitive and straightforward. The actual implementation of the method, if one wants a calculation that runs quickly, involves a little more work. We solve the diffusion equation in Fourier space, where it is diagonal, and back-transform before integrating over the velocity field. With the Neumann boundary conditions, the appropriate Fourier basis is the cosine basis, in which the solution to the diffusion equation has the form

$$\rho (𝐫,t)=\frac{4}{L_xL_y}\underset{𝐤}{}\stackrel{~}{\rho }(𝐤)\mathrm{cos}(k_xx)\mathrm{cos}(k_yy)\mathrm{exp}(k^2t),$$
(8)
where the sum is over all wavevectors $`𝐤=(k_x,k_y)=2\pi (m/L_x,n/L_y)`$ with $`m`$, $`n`$ non-negative integers, and $`\stackrel{~}{\rho }(𝐤)`$ is the discrete cosine transform of $`\rho (𝐫,t=0)`$:

$$\stackrel{~}{\rho }(𝐤)=\{\begin{array}{cc}\frac{1}{4}_0^{L_x}_0^{L_y}\rho (𝐫,0)dxdy\hfill & \text{if }k_x=k_y=0\text{,}\hfill \\ \frac{1}{2}_0^{L_x}_0^{L_y}\rho (𝐫,0)\mathrm{cos}(k_yy)dxdy\hfill & \text{if }k_x=0\text{ and }k_y0\text{,}\hfill \\ \frac{1}{2}_0^{L_x}_0^{L_y}\rho (𝐫,0)\mathrm{cos}(k_xx)dxdy\hfill & \text{if }k_x0\text{ and }k_y=0\text{,}\hfill \\ _0^{L_x}_0^{L_y}\rho (𝐫,0)\mathrm{cos}(k_xx)\mathrm{cos}(k_yy)dxdy\hfill & \text{otherwise.}\hfill \end{array}$$
(9)
The velocity field $`𝐯`$ is then easily calculated from Eqs. (6) and (8) and has components

$`v_x(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{_𝐤k_x\stackrel{~}{\rho }(𝐤)\mathrm{sin}(k_xx)\mathrm{cos}(k_yy)\mathrm{exp}(k^2t)}{_𝐤\stackrel{~}{\rho }(𝐤)\mathrm{cos}(k_xx)\mathrm{cos}(k_yy)\mathrm{exp}(k^2t)}},`$ (10a)
$`v_y(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{_𝐤k_y\stackrel{~}{\rho }(𝐤)\mathrm{cos}(k_xx)\mathrm{sin}(k_yy)\mathrm{exp}(k^2t)}{_𝐤\stackrel{~}{\rho }(𝐤)\mathrm{cos}(k_xx)\mathrm{cos}(k_yy)\mathrm{exp}(k^2t)}}.`$ (10b)
Equations (9) and (10) can be evaluated rapidly using the fast Fourier transform (FFT) and its back-transform respectively, both of which in this case run in time of order $`L_xL_y\mathrm{log}(L_xL_y)`$. We then use the resulting velocity field to integrate Eq. (7), which is a nonlinear Volterra equation of the second kind and can be solved numerically by standard methods . In practice it is the Fourier transform that is the time-consuming step of the calculation and, as we will see, with the aid of the FFT this step can be performed fast enough that the whole calculation runs to completion in a matter of seconds or at most minutes, even for large and detailed maps.

It is straightforward to see that our diffusion cartogram satisfies the fundamental definition, Eq. (1), of a cartogram. In the limit $`t\mathrm{}`$, Eq. (8) is dominated by the $`𝐤=0`$ term and gives

$$\rho (𝐫,\mathrm{})=\frac{4\stackrel{~}{\rho }(0)}{L_xL_y}=\frac{1}{L_xL_y}_0^{L_x}_0^{L_y}\rho (𝐫,0)dxdy=\overline{\rho },$$
(11)
where $`\overline{\rho }`$ is again the mean population density over the area mapped. Furthermore, by definition, the total population within any moving element of area does not change during the diffusion process, and hence, denoting by $`𝐓(𝐫)`$ the final position of a point that starts at position $`𝐫`$, we have $`\rho (𝐫)\mathrm{d}x\mathrm{d}y=\overline{\rho }\mathrm{d}T_x\mathrm{d}T_y`$, and, rearranging, the Jacobian is given by $`(T_x,T_y)/(x,y)=\rho (𝐫)/\overline{\rho }`$, in agreement with Eq. (1).

Below we present a number of applications of our cartograms to the visualization of various types of data. Before we do that however, it is worth giving a brief comparison of our method with the others discussed in Sec. 2. The main advantages of our method over the original method of Tobler , are first speed and second independence of a particular choice of cells. Tobler’s method is obliged to choose directions for the dilation of the cells and this means that the resulting cartogram depends on the coordinate axes used. As we will see, the practical application of our algorithm, like most numerical calculations, relies on a sampling of the population density on a lattice of some kind, but the diffusion process is, in the limit of small lattice spacing, independent of the topology of this lattice and therefore does not show lattice-dependent effects like the method of Tobler. The method of Gusein-Zade and Tikunov  is also free of these effects, but pays for it in run time and complexity.

Conceptually our algorithm is in some respects similar to the cellular automaton method of Dorling . Our description of the diffusion method has been entirely in terms of macroscopic variables and equations, but one could equally look at the method as a microscopic diffusion process in which each individual member of the population performs a Gaussian random walk about the surface of the map. Over time the population will diffuse until it is uniform everywhere within the box enclosing the map, except for statistical fluctuations. The projection of the cartogram is then derived by moving all boundaries on the map in such a way that the net flow passing through them is zero at all times during the diffusion process. This resembles Dorling’s method in the sense that different regions trade their area until a fair distribution is reached.

Our method however has an advantage over Dorling’s in being based on a global, lattice-independent diffusion process. The exchange of area between regions in Dorling’s method occurs only between nearest-neighbor squares along the principle axes of a square lattice and this introduces a strong signature of the lattice topology into the final cartogram.

## 4 Population density function

The description of our method tells, in a sense, only half the story of how to create a cartogram. Before applying this or indeed any method, we need to choose the starting density $`\rho (𝐫)`$ for the map. We can, by defining $`\rho (𝐫)`$ in different ways, control the properties of the resulting cartogram, including crucially the balance between accurate density equalization and readability.

Population density is not strictly a continuous function, since people are themselves discrete and not continuous. To make a continuous function the population must be binned or coarse-grained at some level. All methods of constructing cartograms require one to do this and there is no single accepted standard approach. Part of the art of making a good cartogram lies in shrewd decisions about the definition of the population density.

If we choose a very fine level of coarse-graining for the population density, then the high populations in centers such as cities will require substantial local distortions of the map in order to equalize the density. A coarser population density will cause less distortion, resulting in a map with features that are easier to recognize, but will give a less accurate impression of the true population distribution. The most common choice made by others has been to coarse-grain the population at the level of the (usually political) regions of interest. For example, if one were interested in the United States, one might take the population of each state and distribute it uniformly over the area occupied by that state. This method can be used also with our cartogram algorithm and we give one example below. But we are not obliged to use it, and in some cases it may be undesirable, since binning at the level of states erases any details of population distribution below the state level. In our work we have mostly adopted instead a spatially uniform Gaussian blur as our coarse-graining function. By varying the width $`\sigma `$ of the blurring function we can tune our cartogram between accuracy and readability.

Ultimately the choice of population density function is up to the user of the method, who must decide what particular features are most desirable in his or her application. One advantage of our diffusion-based algorithm is that it is entirely agnostic about this choice; the process of computing the cartogram is completely decoupled from the calculation of the population density and hence is not slanted in favor of one choice or another.

## 5 Applications

We give three examples of the use of our cartograms, focusing on the United States and using population data from the 2000 US Census.

First, we examine the results of the US presidential election of 2000. Figure 1a shows the popular vote by state in this election—the simple fraction of voters voting each way—with state majorities for the two major candidates represented as shades of red (George Bush, Republican) and blue (Al Gore, Democrat). The briefest appraisal by eye immediately reveals that the Republicans dominate much more than a half of the country. This, however, is a misleading conclusion because much of the Republicans’ dominance comes from their success in the large but relatively unpopulated states in the center of the map, while the Democrats carry the more populous areas in the northeast and on the west coast. Clearly then, this is a poor visual representation of the results, in the sense that it is hard to tell which party got more votes by looking at the map.

A much better representation is given by Figure 1b, in which the same color code is used on a population cartogram constructed using the methods described in this paper. This cartogram is based on a density function with a moderately broad Gaussian blur, as described above, yielding a map whose features are distorted relatively little on short scales; the individual states are still easily recognizable while being scaled close to the size appropriate to their populations. To a good approximation the amounts of red and blue in the figure now correspond to the true balance of the popular vote, and, as is clear to the eye, this vote was very close between the two parties, in fact being won not by the Republican candidate but by the Democrat. In Fig. 1c, we show a further cartogram constructed using less coarse-graining on the population density, resulting in a map that more perfectly represents states’ populations, but which also has more distortion on short length-scales, making some regions hard to recognize. For example, the small but densely populated Long Island now expands (quite correctly) to a size greater than the entire state of Wyoming. The user concerned both with readability and accurate portrayal of the data would probably choose a map similar to Fig. 1b in this case.

Ultimately, the presidency is decided not by the popular vote, but by the electoral college. Under the US system, each state contributes a certain number of electors to the electoral college, who vote according to the majority in their state. The candidate receiving a majority of the votes in the electoral college wins the election. The appropriate visualization for a vote of this kind is one in which the sizes of the states are scaled in proportion to their numbers of electors. This then is an example in which a coarse-graining according to political boundaries (state boundaries in this case) makes good sense. We show a cartogram calculated in this way in Fig. 1d. The allocation of electors to states roughly follows population levels, but contains a deliberate bias in favor of less populous states, and as a result some of these states appear larger in Fig. 1d than in 1c—Wyoming, Montana, and the Dakotas are good examples. Since most of these are majority-Republican states, we can now understand how the Republican candidate came to win the election in spite of losing the popular vote.

For our second example, we look at a case in which a very fine level of coarse-graining is needed to understand the data fully. We study the distribution of cases of lung cancer among the male population in the state of New York. (A similar study using a different technique and for a smaller area was carried out by Merrill .) In Fig. 2a we show a scatter plot of lung cancer cases between 1993 and 1997, as recorded by the New York State Department of Health. This map is of precisely the kind discussed in the introduction: it shows clearly how many cases there are, and that they are geographically concentrated in the areas that have high populations. However, it is impossible to tell whether there is a statistically higher per capita incidence of lung cancer in one area or another, because any such variation is masked by the highly non-uniform population density.

In Fig. 2b, we show the same data on a population cartogram with moderate coarse-graining of the initial population density. Although the map is visibly distorted, little difference is visible in the distribution of cancer cases. In Fig. 2c, on the other hand, we use a very fine-grained population density, creating a cartogram with better population equalization and significantly greater distortion. Now the virtue of this representation becomes strikingly clear. As the figure shows, when we use a projection that truly equalizes the population density over the map, there is no longer any significant variation in the distribution of cases over the state—the dots have about the same density everywhere. The shape of the map in (c) does not much resemble the shape of the original any more, but this is the price we pay for equalizing the population almost perfectly.

Our method of generating cartograms is fast, an important consideration for interactive use. As discussed in Sec. 3, the bulk of the work involved in creating the maps is in the Fourier transforms, which can be computed rapidly using FFTs. Fig. 2c for example was produced in less than 100 seconds on a standard desktop computer, including the time to read in the census data, to perform the Gaussian convolution, to solve the diffusion equation, and to plot the figure. Previous techniques are either significantly slower (Kocmoud  reports 16 hours for a US state cartogram using his constraint-based approach) or are obliged to use simplified maps to reduce the computational load .

The example cartograms given so far have all been based, more or less, on human population density, which is certainly the most common type of cartogram. Other types however are also possible and for our third example we study one such. Anyone who reads or watches the news in the United States (and similar observations probably apply in other countries as well) will have noticed that the geographical distribution of news stories is skewed. Even allowing for population, a few cities, notably New York and Washington, DC, get a surprisingly large fraction of the attention while other places get little. Apparently some locations loom larger in our mental map of the nation than others, at least as presented by the major media. We can turn this qualitative idea of a mental map into a real map using our cartogram method.

We have taken about $`\mathrm{72\hspace{0.17em}000}`$ newswire stories from November 1994 to April 1998, and extracted from each the “dateline,” a line at the head of the story that gives the date and the location that is the main focus of the story. Binning these locations by state, we then produce a “mindshare map” in which the sizes of the US states represent the fraction of stories concerning that state over the time interval in question. The result is shown in Fig. 3.

The stories are highly unevenly distributed. New York City alone contributes $`\mathrm{20\hspace{0.17em}000}`$ stories to the corpus—largely because of the preponderance of wire stories about the financial markets—and Washington, DC another $`\mathrm{10\hspace{0.17em}000}`$, largely political stories. We chose to bin by state to avoid large distortions around the cities that are the focus of most news stories. We made one exception however: since New York City had far more hits than any other location including the rest of the state of New York (which had around 1000), we decided to split New York State into two regions, one for the greater New York City area and another for the rest of the state.

The cartogram is a dramatic depiction of the distribution of US news stories. The map is highly distorted because the patterns of reporting show such extreme variation. Washington, DC, for instance, which is normally invisible on a map of this scale, becomes the second largest “state” in the union. (The District of Columbia is not, technically, a state.) Members of the public, including the US public, frequently overestimate the size of the northeastern part of the United States by comparison with the middle and western states, and this map may give us a clue as to why. Perhaps people’s mental image of the United States is not really an inaccurate one; it is simply based on things—such as “mindshare”—that people find more relevant than geographical area.

## 6 Conclusions

In this paper we have presented a new general method for constructing density-equalizing projections or cartograms, which provide an invaluable tool for the presentation of geographic data. Our method is simpler than many earlier methods, allowing for rapid calculations, while generating accurate and readable maps. The method allows its users to choose their own balance between good density equalization and low distortion of map regions, making it flexible enough for a wide variety of applications. We have given comparisons between our method and previous ones, and presented a number of examples of the use of our cartograms in the representation of human data, including election results and incidence data for cancer.

Implementation of our method in GIS software packages should be straightforward, and we hope that in this or other forms it will prove a valuable tool for researchers in a wide variety of disciplines.

## Acknowledgments

The authors would like to thank the staff of the University of Michigan’s Numeric and Spatial Data Services for their help with the geographic data, and Dragomir Radev for useful discussions about the geographic distribution of news messages. This work was funded in part by the National Science Foundation under grant number DMS–0234188 and by the James S. McDonnell Foundation.


text:
# Measuring polynomial functions of states

## I Introduction

Much of the practice of quantum information theory revolves around calculations of various information measures. These are functions of the quantum state $`\rho `$ which quantify properties of interest, such as purity, entanglement, and distinguishability. An important subset of these are polynomial functions. Some polynomial quantities of potential interest are unitary invariants such as Kempe’s invariant and the 3-tangle Kempe99 ; Coffman00 ; Carteret03a , and measures of entanglement such as the $`Q`$ measure of Meyer and Wallach Meyer01 ; Brennen03 . Some other quantities, which are not polynomials themselves, can either be approximated by polynomials (such as the quadratic approximation of the von Neumann entropy) or are simple functions of such polynomials, such as the concurrence and the negativity of the partial transpose Wootters98 ; Peres96 ; Carteret03b . (I state what is meant by a polynomial function of $`\rho `$ more precisely below.)

A difficulty in estimating polynomial functions of the state is the linearity of quantum mechanics. Given only one copy of a system in a particular quantum state $`\rho `$, any measurement can only depend linearly on $`\rho `$. If it is possible to repeatedly prepare systems in the state $`\rho `$, we can overcome this obstacle by performing quantum state tomography Vogel89 : by estimating the expectation values in the state $`\rho `$ of a set of observables, one can construct a classical description of the state; given that description, obviously any function of $`\rho `$ can be calculated. However, many measurements are needed, of many different observables, since the number of parameters to be estimated grows like $`d^2`$ where $`d`$ is the dimension of the system’s Hilbert space.

Recently there has been considerable interest in the problem of determining polynomial functions of quantum states (or simple functions of such polynomials); and a number of ingenious procedures have been put forward for measuring particular functions Sjoqvist00 ; Keyl01 ; Horodecki01a ; Horodecki01b ; Horodecki01c ; Fiurasek01 ; Ekert02 ; Filip02 ; Alves03 . These procedures try, by various clever techniques, to avoid estimating all the parameters needed for full state tomography; rather, they either estimate the polynomial function directly, or estimate a smaller set of parameters from which the polynomial can be calculated. They do this by performing measurements on several copies of $`\rho `$ at once.

In this paper I demonstrate that for any $`m`$th-degree polynomial $`f`$ in the elements of the density matrix $`\rho `$ (on a Hilbert space $``$), there are observables $`\widehat{O}_f`$ and $`\widehat{O}_f^{}`$ on $`^m`$ such that the expectation value of $`\widehat{O}_f+i\widehat{O}_f^{}`$ on $`m`$ copies of $`\rho `$ equals $`f`$:

$$f=\mathrm{Tr}\{(\widehat{O}_f+i\widehat{O}_f^{})\rho ^m\}\widehat{O}_f+i\widehat{O}_f^{}.$$
(1)
## II Polynomial functions

What do I mean by a polynomial function on the elements of $`\rho `$? Let us suppose that $`\{|i\}`$ is a basis for the Hilbert space $``$, where $`0i<d`$. Then we can write $`\rho `$ in terms of that basis.

$$\rho =\underset{i,j=0}{\overset{d1}{}}\rho _{ij}|ij|.$$
(2)
A polynomial function $`f(\rho )`$ of degree $`m`$ can be written

$$f(\rho )=\underset{i_1,j_1,\mathrm{},i_m,j_m}{}c_{i_1j_1\mathrm{}i_mj_m}\rho _{i_1j_1}\rho _{i_2j_2}\mathrm{}\rho _{i_mj_m},$$
(3)
where the $`\{c_{i_1j_1\mathrm{}i_mj_m}\}`$ are arbitrary complex constants. There could also be terms of lower-order than $`m`$; but we can exploit the fact that $`\mathrm{Tr}\rho =_i\rho _{ii}=1`$ to write them all as sums of terms of order $`m`$, so the form (3) is general. This is easy to see: consider a term of order $`k<m`$, proportional to $`\rho _{i_1j_1}\mathrm{}\rho _{i_kj_k}`$. We can rewrite this as a sum of terms of order $`m`$ by multiplying it by $`mk`$ factors of $`\mathrm{Tr}\rho `$:

$$\rho _{i_1j_1}\mathrm{}\rho _{i_kj_k}=\underset{i_{k+1},\mathrm{},i_m}{}\rho _{i_1j_1}\mathrm{}\rho _{i_kj_k}\rho _{i_{k+1}i_{k+1}}\mathrm{}\rho _{i_mi_m}.$$
(4)
General polynomials of form (3) will not be invariant under a change of basis (though some polynomials of interest are invariant). However, since a basis change is a linear transformation, such a transformation takes $`m`$th-degree polynomials to $`m`$th-degree polynomials.

Suppose we have $`m`$ copies of the system in state $`\rho `$; these copies have the joint state

$$\rho ^m=\underset{i_1,j_1,\mathrm{},i_m,j_m}{}\rho _{i_1j_1}\mathrm{}\rho _{i_mj_m}|i_1j_1|\mathrm{}|i_mj_m|.$$
(5)
Consider one term of the polynomial function $`f(\rho )`$ given by (3), proportional to $`\rho _{i_1j_1}\mathrm{}\rho _{i_mj_m}`$. We can write down an operator

$$\widehat{A}_{i_1j_1\mathrm{}i_mj_m}=|j_1i_1|\mathrm{}|j_mi_m|,$$
(6)
such that

$$\mathrm{Tr}\{\widehat{A}_{i_1j_1\mathrm{}i_mj_m}\rho ^m\}=\rho _{i_1j_1}\mathrm{}\rho _{i_mj_m},$$
(7)
as we see by substituting (6) and (5) in (7). Since the full polynomial $`f(\rho )`$ is a linear combination of such terms, and the trace is a linear operation, we can find an operator $`\widehat{A}_f`$

$$\widehat{A}_f=\underset{i_1,j_1,\mathrm{},i_d,j_d}{}c_{i_1j_1\mathrm{}i_mj_m}\widehat{A}_{i_1j_1\mathrm{}i_mj_m},$$
(8)
such that

$$f(\rho )=\widehat{A}_f=\mathrm{Tr}\{\widehat{A}_f\rho ^m\}.$$
(9)
While this operator obviously exists, it will in general not be an observable—that will only be the case if $`\widehat{A}_f`$ is also Hermitian, $`\widehat{A}_f=\widehat{A}_f^{}`$. However, we can trivially find a pair of operators which are observables,

$`\widehat{O}_f`$ $`=`$ $`(\widehat{A}_f+\widehat{A}_f^{})/2,`$
$`\widehat{O}_f^{}`$ $`=`$ $`i(\widehat{A}_f\widehat{A}_f^{})/2,`$ (10)
such that $`\widehat{A}_f=\widehat{O}_f+i\widehat{O}_f^{}`$.

These operators will not be unique, in general. If $`\widehat{A}_f`$ is an operator as defined above, then the operator obtained by permuting the order of the $`m`$ systems will also give the polynomial. So will any linear combination of such permutations, so long as the coefficients of the linear combination sum to 1. This freedom reflects the permutation symmetry of the state $`\rho ^m`$

## III Measuring polynomial functions

If we have the ability to measure arbitrary observables, then we are done; we repeatedly prepare $`\rho ^m`$ and measure $`\widehat{O}_f`$ and $`\widehat{O}_f^{}`$ until we have a good estimate of their expectation values. If we do not have the ability to measure arbitrary observables, the fact that an observable exists whose expectation value yields the polynomial function is not very useful in itself. However, even if we cannot measure arbitrary observables, we can design a circuit which will simulate such a measurement to any desired degree of accuracy.

The proof of this is trivial. Let the observable $`\widehat{O}_f`$ have a basis of eigenstates $`|\varphi _j`$ with corresponding eigenvalues $`o_j`$ (which are not necessarily distinct). Let $`|j`$ be a standard basis which we are capable of measuring. (For instance, in a system of q-bits this could be the computational basis.) Then we can define a unitary transformation

$$\widehat{U}_f=\underset{j}{}|j\varphi _j|,$$
(11)
which rotates the observable $`\widehat{O}_f`$ to be diagonal in the standard basis:

$$\widehat{O}_f=\underset{j}{}o_j|\varphi _j\varphi _j|,\widehat{U}_f\widehat{O}_f\widehat{U}_f^{}=\underset{j}{}o_j|jj|.$$
(12)
Clearly this new observable $`\widehat{U}_f\widehat{O}_f\widehat{U}_f^{}`$ is straightforward to measure: we just measure the state in the standard basis and give a value $`o_j`$ to outcome $`j`$. This gives us an expression for the expectation value of $`\widehat{O}_f`$,

$`\mathrm{Tr}\{\widehat{O}_f\rho ^m\}`$ $`=`$ $`\mathrm{Tr}\{(\widehat{U}_f\widehat{O}_f\widehat{U}_f^{})(\widehat{U}_f\rho ^m\widehat{U}_f^{})\}`$ (13)
$`=`$ $`{\displaystyle \underset{j}{}}o_jj|(\widehat{U}_f\rho ^m\widehat{U}_f^{})|j.`$
All we need then is to be able to find a circuit which approximates the unitary transformation $`\widehat{U}_f`$. It is well-known that a circuit exists to approximate any unitary transformation to any desired degree of accuracy Divincenzo95 . Our procedure now is as follows: repeatedly prepare $`m`$ copies of the state $`\rho `$, and carry out the circuit which approximates $`\widehat{U}_f`$. Measure in the standard basis, and weight each outcome $`j`$ with the eigenvalue $`o_j`$. Do the same for the observable $`\widehat{O}_f^{}`$. We see, therefore, that in principle it is always possible to find a circuit which enables one to measure any polynomial function of $`\rho `$.

## IV Finding circuits

This proof of principle is very far from being a proof that such a measurement is practical. While such a circuit always exists, there is no guarantee that it will achieve the desired unitary transformation efficiently. Indeed, most unitary transformations cannot be approximated efficiently Knill95 . A protocol of this nature would only be worthwhile if it were substantially more efficient than performing quantum state tomography on $`\rho `$ and then calculating the function directly. This requires comparing the cost of doing more measurements (in the case of state tomography) to the cost of preparing a more complicated state ($`\rho ^m`$ instead of $`\rho `$) and performing some (possibly large) number of extra gates.

Even if it is possible, one must in most cases diagonalize a pair of complicated observables $`\widehat{O}_f`$ and $`\widehat{O}_f^{}`$ and find the corresponding eigenvalues and eigenvectors. This may be a nontrivial task; if the system is of dimension $`d`$, then $`\widehat{O}_f`$ is a $`d^m\times d^m`$ Hermitian matrix. However, for some applications this might be acceptable; the eigenvectors and eigenvalues could be calculated “off-line,” and then used to design a circuit which could be applied to many different states $`\rho `$. There is no known efficient algorithm for finding the simplest circuit that produces a given unitary transformation. Indeed, that is almost certain a computationally intractable problem in itself. But in some cases, it may be possible to find a simple circuit to do the job. Polynomials with a great deal of symmetry may simplify the circuits considerably.

Let us consider a simple example. The first, and undoubtedly most immediately useful, is calculating the trace of $`\rho ^2`$:

$$f=\mathrm{Tr}\{\rho ^2\}=\underset{i,j}{}\rho _{ij}\rho _{ji}=\underset{i,j}{}|\rho _{ij}|^2.$$
(14)
This quantity gives a useful measure of the purity of a state, with $`\mathrm{Tr}\{\rho ^2\}=1`$ for a pure state and $`\mathrm{Tr}\{\rho ^2\}=1/d`$ for the maximally mixed state; this is easier to calculate than the von Neumann entropy. For a bipartite pure state $`|\mathrm{\Psi }`$, $`\mathrm{Tr}\{\rho _A^2\}`$ gives a measure of the two subsystems’ entanglement, where $`\rho _A=\mathrm{Tr}_B\{|\mathrm{\Psi }\mathrm{\Psi }|\}`$ is the reduced density matrix for subsystem A. From (14) it is obvious that the correct observable is the swap operator:

$$\widehat{O}_f=\underset{i,j}{}|ij||ji|.$$
(15)
This is a special case of the widely known fact that

$$\mathrm{Tr}\{\rho ^m\}=\mathrm{Tr}\{\widehat{S}\rho ^m\}$$
(16)
where $`\widehat{S}`$ is the cyclic shift operator: $`\widehat{S}|\psi _1\mathrm{}|\psi _m=|\psi _m|\psi _1\mathrm{}|\psi _{m1}`$. In the case $`m=2`$, this is the pairwise swap, and is Hermitian as well as unitary, so that it is itself an observable. It has eigenvalues $`\pm 1`$, as any operator which is both Hermitian and unitary must; and if our system is a single q-bit, a good eigenbasis is the usual Bell basis:

$`|\varphi _0`$ $`=`$ $`(|00+|11)/\sqrt{2},o_0=1,`$
$`|\varphi _1`$ $`=`$ $`(|00|11)/\sqrt{2},o_1=1,`$
$`|\varphi _2`$ $`=`$ $`(|01+|10)/\sqrt{2},o_2=1,`$
$`|\varphi _3`$ $`=`$ $`(|01|10)/\sqrt{2},o_3=1.`$ (17)
Circuits for measuring the Bell states are widely known; Fig. 1 presents an example of such a circuit, consisting of a controlled-NOT gate, a Hadamard gate, and two single-bit meaurements in the computational basis. In many cases it would be much easier to measure such an observable than to do full state tomography even on a single bit.

## V Relationship to other approaches

Other circuits have been presented in the literature for calculating other polynomial functions of the state Horodecki01b ; Ekert02 ; Filip02 ; Alves03 ; Paz03 . It is not difficult to see that these circuits give a method for estimating the expectation value of a suitable observable $`\widehat{O}_f`$, as described in this paper.

Let us look at this a bit more closely. Consider the circuit in Fig. 2, which is a version of the circuit presented in Ekert02 . In this circuit, a set of $`m`$ systems is prepared in the state $`\rho ^m`$; an extra “control bit” is prepared in the state $`|0`$. This control bit undergoes a Hadamard gate; it then serves as the control for a controlled-$`\widehat{U}`$ gate (or circuit) with the $`m`$ systems, followed by another Hadamard. Finally, the observable $`\widehat{Z}=|00||11|`$ is measured for the control bit. It is straightforward to show that the expectation value of this measurement is $`\widehat{Z}=(1/2)\mathrm{Tr}\{(\widehat{U}+\widehat{U}^{})\rho ^m\}`$. If $`\widehat{U}=\widehat{S}`$, for instance, then this would be equal to $`\mathrm{Tr}\{\rho ^m\}`$. The operator $`(1/2)(\widehat{U}+\widehat{U}^{})`$ is Hermitian and therefore an observable.

Why does this work? The values of the measurements can only be $`1`$ and $`1`$; for most observables these will not be equal to any eigenvalue of $`(1/2)(\widehat{U}+\widehat{U}^{})`$. In spite of this, the expectation values are exactly the same as if the observable $`(1/2)(\widehat{U}+\widehat{U}^{})`$ had been measured. To see what is going on, let us suppose for the moment that $`\widehat{U}`$ has eigenstates $`|\varphi _j`$ with eigenvalues $`\mathrm{exp}(i\theta _j)`$. We can always express the state of the $`m`$ systems in terms of these eigenstates,

$$\rho ^m=\underset{j,j^{}}{}R_{jj^{}}|\varphi _j\varphi _j^{}|.$$
(18)
Suppose that we start the control bit and the $`m`$ systems in the joint state $`|\mathrm{\Psi }=|0|\varphi _j`$, and carry out the circuit in Fig. 2. The resulting state will be

$$|\mathrm{\Psi }^{}=\mathrm{exp}(i\theta _j/2)\left[\mathrm{cos}(\theta _j/2)|0i\mathrm{sin}(\theta _j/2)|1\right]|\varphi _j\mathrm{exp}(i\theta _j/2)|\gamma _j|\varphi _j.$$
(19)
The expectation $`\widehat{Z}`$ in the control-bit state $`|\gamma _j`$ is $`\mathrm{cos}^2(\theta _j/2)\mathrm{sin}^2(\theta _j/2)=\mathrm{cos}(\theta _j)`$ which is the $`j`$th eigenvalue of the operator $`(1/2)(\widehat{U}+\widehat{U}^{})`$. A similar result is obtained for each eigenstate $`|\varphi _j`$. This circuit is therefore equivalent to the unitary transformation

$$\widehat{E}=\underset{j}{}\mathrm{e}^{i\theta _j/2}\mathrm{exp}(i\theta _j\widehat{X}/2)|\varphi _j\varphi _j|.$$
(20)
So a different way of understanding this circuit sees the control bit as a target bit: if the $`m`$ systems are in the $`j`$th eigenstate of $`\widehat{U}`$, the target bit is rotated about the $`x`$ axis from $`|0`$ to a new state $`|\gamma _j`$ whose $`\widehat{Z}`$ expectation exactly matches the $`j`$th eigenvalue of the observable $`(1/2)(\widehat{U}+\widehat{U}^{})`$. If we now start the $`m`$ systems in the state $`\rho ^n`$, the different eigenstates of $`\widehat{U}`$ will contribute with their appropriate weights $`R_{jj}`$:

$`\mathrm{Tr}\{(\widehat{Z}\widehat{I}^m)\widehat{E}(|00|\rho ^m)\widehat{E}^{}\}`$ $`=`$ $`{\displaystyle \underset{jj^{}}{}}R_{jj^{}}\gamma _j^{}|\widehat{Z}|\gamma _j\varphi _j^{}|\varphi _j`$ (21)
$`=`$ $`{\displaystyle \underset{j}{}}R_{jj}\gamma _j|\widehat{Z}|\gamma _j`$
$`=`$ $`{\displaystyle \underset{j}{}}R_{jj}\varphi _j|(1/2)(\widehat{U}+\widehat{U}^{})|\varphi _j`$
$`=`$ $`{\displaystyle \underset{jj^{}}{}}R_{jj^{}}\varphi _j^{}|(1/2)(\widehat{U}+\widehat{U}^{})|\varphi _j`$
$`=`$ $`\mathrm{Tr}\{(1/2)(\widehat{U}+\widehat{U}^{})\rho ^m\}.`$
Clearly, for any observable $`\widehat{O}`$ whose eigenvalues $`o_j`$ lie in the range $`1o_j1`$, we can find a unitary $`\widehat{U}`$ such that $`\widehat{O}=(1/2)(\widehat{U}+\widehat{U}^{})`$. For example, $`\widehat{U}=\widehat{O}+i\sqrt{\widehat{I}\widehat{O}^2}`$ is such an operator. If the eigenvalues of $`\widehat{O}`$ lie outside this range, we can construct a new observable $`\widehat{O}/c`$ where $`c>0`$ is larger than the absolute value of the largest eigenvalue of $`\widehat{O}`$. We can then find a unitary $`\widehat{U}`$ such that $`(1/2)(\widehat{U}+\widehat{U}^{})=\widehat{O}/c`$, and estimate the expectation value of $`\widehat{O}/c`$ by carrying out the circuit in Fig. 2. This is an alternative and in some cases much more efficient way of estimating a polynomial function of $`\rho `$; in particular, it only requires measurements of a single bit. Paz and Roncaglia Paz03  have constructed a programmable circuit of this type to estimate the expectation value of any observable. This would obviously be sufficient to estimate any polynomial function of the state $`\rho `$.

## VI Conclusions

The ability to measure arbitrary $`m`$-system observables and determine their expectation values allows the estimation of arbitrary $`m`$th-degree polynomials of the state, without having to perform quantum state tomography. This is an enormous simplification, in principle; in some cases, it may give a considerable simplification in practice. Whatever the practicalities, this equivalence is a useful fact, which may make possible improved protocols for estimating entanglement measures, polynomial approximations to the entropy, or other functions of the state.

###### Acknowledgements.

I would like to acknowledge financial support from the Martin A. and Helen Chooljian Membership in Natural Sciences, and DOE Grant No. DE-FG02-90ER40542. I appreciate helpful feedback on these ideas from Hilary Carteret and Rüdiger Schack.


text:
# Boyer-Finley equation and systems of hydrodynamic type

Ferapontov E.V.,

Department of Mathematical Sciences

Loughborough University

Loughborough, Leicestershire LE11 3TU

United Kingdom

e-mail: E.V.Ferapontov@lboro.ac.uk
Korotkin D.A., Shramchenko V.A.
Department of Mathematics and Statistics

Concordia University

7141, rue Sherbrooke Ouest

Montréal, Québec H4B 1R6

Canada

e-mails: korotkin@mathstat.concordia.ca
vasilisa@alcor.concordia.ca

Abstract. We reduce Boyer-Finley equation to a family of compatible systems of hydrodynamic type, with characteristic speeds expressed in terms of spaces of rational functions. The systems of hydrodynamic type are then solved by the generalized hodograph method, providing solutions of the Boyer-Finley equation including functional parameters.

In this paper we construct solutions of the dispersionless non-linear PDE – the Boyer-Finley equation (self-dual Einstein equation with a Killing vector),

$$U_{xy}=(e^U)_{tt},$$
(1)
via reduction to a family of compatible systems of hydrodynamic type.

This equation was actively studied during last twenty years by many authors; we just mention works . So far the most general scheme of the construction of its solutions was developed in . In these works solutions of the Boyer-Finley equation were derived by averaging an appropriate two-point Baker-Akhiezer function in genus zero which corresponds to the two-dimensional Toda lattice equations, the underlying Riemann surface deforming according to the Whitham equations. Some particular solutions of the Boyer-Finley equation were constructed in ; their relationship to the solutions of  remains unclear.

The goal of the present paper is to give an alternative scheme of solving the Boyer-Finley equation. Namely, we consider reductions of this equation to multi-component systems of hydrodynamic type in the spirit of , see also ; the equations for charactersitic speeds of these systems are solved in terms of rational branched coverings. The systems of hydrodynamic type are then solved by the generalized hodograph method .

Now we describe a way to find solutions of the Boyer-Finley equation. First we ignore the condition of reality of the function $`U`$ and construct complex-valued solutions of equation (1); then we formulate restrictions on the parameters which guarantee the reality of $`U.`$

Let us assume $`U`$ to be a function of $`L`$ variables $`\lambda _1,\mathrm{},\lambda _L`$, where $`\lambda _m(x,y,t)`$ (“Riemann invariants”) satisfy a pair of systems of hydrodynamic type

$$_x\lambda _m=V_m(\{\lambda _k\})_t\lambda _m,_y\lambda _m=W_m(\{\lambda _k\})_t\lambda _m.$$
(2)
A direct substitution of the function $`U(\lambda _1,\mathrm{},\lambda _L)`$ into the Boyer-Finley equation implies the algebraic relation among the functions $`U,V_m`$ and $`W_m`$,

$$V_mW_m=e^U,$$
(3)
along with the following differential equations:

$$_m_nU(V_mW_n+V_nW_m)=2_n_m(e^U),mn,$$
(4)
$$\frac{_mV_n}{V_mV_n}=\frac{_mW_n}{W_mW_n},mn,$$
(5)
where $`_m/\lambda _m`$.

Relations (3) allow one to parameterize the functions $`V_m`$ and $`W_m`$ by a new set of variables $`\phi _m(\{\lambda _k\})`$ as follows:

$$V_m=\mathrm{exp}\left\{2i\phi _m+\frac{U}{2}\right\},W_m=\mathrm{exp}\left\{2i\phi _m+\frac{U}{2}\right\}.$$
(6)
In terms of these variables, equations (4) and (5) take the form

$$_m_nU=\frac{_mU_nU}{2\mathrm{sin}^2(\phi _m\phi _n)},_m\phi _n=\frac{1}{4}\mathrm{cot}(\phi _n\phi _m)_mU,$$
(7)
where $`mn`$.

A class of solutions of this system is related to the space of rational functions in the following way. Consider a rational function

$$R(\mu )=\mu +\underset{k=1}{\overset{N1}{}}\frac{a_k}{\mu b_k},\mu ^1.$$
(8)
In application to Benney’s hierarchy, functions of this form were first introduced in . The equation

$$\lambda =R(\mu )$$
(9)
defines an $`N`$sheeted covering $``$ of the $`\lambda `$sphere. A point $`P`$ is a pair of complex numbers, $`P=(\lambda ,\mu )`$. We consider the generic case when the function $`R(\mu )`$ has $`2N2`$ non-coinciding finite critical points, i.e., the equation

$$R^{}(\mu )=0$$

has $`2N2`$ distinct roots $`\mu _1,\mathrm{},\mu _{2N2}.`$ The corresponding critical values,

$$\lambda _n=R(\mu _n),n=1,\mathrm{},2N2,$$
(10)
are projections onto the $`\lambda `$sphere of the branch points of the covering $``$ (we denote branch points by $`P_n=(\lambda _n,\mu _n)`$; all of them are simple as a corollary of non-coincidence of $`\mu _n`$ for different $`n`$). An additional condition we impose on the function $`R(\mu )`$ is that all $`\lambda _n`$ are different. Now, observe that the number of parameters of the rational function (8) is equal to the number of branch points; therefore, we can take $`\lambda _1,\mathrm{},\lambda _{2N2}`$ as local coordinates on the space of rational functions. It was shown in  that the critical points $`\{\mu _m\}`$ of the rational function $`R(\mu )`$ depend on $`\{\lambda _n\}`$ in the following way:

$`{\displaystyle \frac{\mu _m}{\lambda _n}}={\displaystyle \frac{\beta _n}{\mu _n\mu _m}};{\displaystyle \frac{\beta _m}{\lambda _n}}={\displaystyle \frac{2\beta _n\beta _m}{(\mu _n\mu _m)^2}},mn.`$ (11)
These equations appeared also in  in the theory of hydrodynamic reductions of Benney’s moment equations. The inverse function $`\mu (P)=R^1(P)`$ is defined on the covering $``$ . As a function of $`\{\lambda _n\}`$, it satisfies the system of differential equations :

$$\frac{\mu }{\lambda _n}=\frac{\beta _n}{\mu _n\mu }.$$
(12)
Let us now choose two points $`Q_1`$ and $`Q_2`$ on the covering $``$ such that their projections $`\lambda (Q_1)`$ and $`\lambda (Q_2)`$ onto the $`\lambda `$sphere do not depend on $`\{\lambda _n\}`$. Then, consider the following function:

$$\gamma (P)=\frac{1}{2\pi i}\mathrm{log}\frac{\mu (P)\mu (Q_1)}{\mu (P)\mu (Q_2)},$$
(13)
which maps $``$ onto a cylinder. We shall be interested in the images $`\{\gamma _m\}`$ of branch points under this map as functions of $`\{\lambda _n\}.`$ In the sequel we denote $`\mu (Q_1)`$ by $`\kappa _1`$ and $`\mu (Q_2)`$ by $`\kappa _2`$. According to (12), they satisfy the equations

$$\frac{\kappa _j}{\lambda _n}=\frac{\beta _n}{\mu _n\kappa _j},j=1,2.$$
(14)
From the expression (13) for $`\gamma (P_m)`$ we have

$$\gamma _m=\frac{1}{2\pi i}\mathrm{log}\frac{\mu _m\kappa _1}{\mu _m\kappa _2}.$$
(15)
Differentiation of this relation with respect to $`\lambda _n`$ using (11), (13) and (14) gives

$`{\displaystyle \frac{\gamma _m}{\lambda _n}}={\displaystyle \frac{1}{2\pi i}}{\displaystyle \frac{\beta _n}{\mu _n\mu _m}}\left[{\displaystyle \frac{1}{\mu _m\kappa _1}}{\displaystyle \frac{1}{\mu _m\kappa _2}}\right].`$ (16)
If we now express $`\mu _n`$ and $`\mu _m`$ in terms of $`\gamma _n,\gamma _m`$ from (15) and set

$$\alpha _n=\frac{\beta _n}{4\pi ^2}\left[\frac{1}{\mu _m\kappa _1}\frac{1}{\mu _m\kappa _2}\right]^2,$$
(17)
we arrive at the following system of differential equation for the functions $`\gamma _m(\{\lambda _n\})`$:

$$\frac{\gamma _m}{\lambda _n}=\pi \alpha _n\left(\mathrm{cot}\pi (\gamma _m\gamma _n)+\mathrm{cot}\pi \gamma _n\right),mn.$$
(18)
Similarly, the functions $`\alpha _n`$ satisfy the following equations:

$$\frac{\alpha _n}{\lambda _m}=2\pi ^2\frac{\alpha _n\alpha _m}{\mathrm{sin}^2\pi (\gamma _n\gamma _m)},mn.$$
(19)
It turns out that a simple transformation allows one to construct solutions of system (7) from the set of functions $`\gamma _m`$ and $`\alpha _m`$. Namely, the system of equations (18), (19) coincides with system (7) if rewritten in terms of the new variables $`U`$ and $`\phi _n`$ such that

$$\frac{U}{\lambda _m}=4\pi ^2\alpha _m$$
(20)
and

$$\phi _n=\pi (\gamma _n+\psi ),$$
(21)
where

$$\frac{\psi }{\lambda _m}=\pi \alpha _m\mathrm{cot}\pi \gamma _m.$$
(22)
The existence of functions $`U`$ and $`\psi `$ is provided by the compatibility conditions,

$$\frac{}{\lambda _n}\alpha _m=\frac{}{\lambda _m}\alpha _n$$
(23)
and

$$\frac{}{\lambda _n}\left(\alpha _m\mathrm{cot}\pi \gamma _m\right)=\frac{}{\lambda _m}\left(\alpha _n\mathrm{cot}\pi \gamma _n\right),$$
(24)
which follow from (19) and (18).

Ultimately, formulae (6) determine $`\{V_m\}`$ and $`\{W_m\}`$ as functions of $`\{\lambda _k\}`$. In order to obtain a solution of the Boyer-Finley equation (1) we need $`U`$ as an explicit function of $`x`$, $`y`$ and $`t`$, that is, we need to solve the system of hydrodynamic type (2). The tool which is usually used for this purpose is the generalized hodograph method . Instead of solving (2), we find a smooth solution $`\mathrm{\Lambda }(x,y,t)=(\lambda _1,\mathrm{},\lambda _L)`$ of the following system

$$\mathrm{\Phi }_m(\mathrm{\Lambda })=t+V_m(\mathrm{\Lambda })x+W_m(\mathrm{\Lambda })y,$$
(25)
where the functions $`\{\mathrm{\Phi }_m\}`$ satisfy the linear system

$$\frac{_m\mathrm{\Phi }_n}{\mathrm{\Phi }_m\mathrm{\Phi }_n}=\frac{_mV_n}{V_mV_n}=\frac{_mW_n}{W_mW_n},mn.$$
(26)
To see that an implicit solution (25) for $`\{\lambda _m(x,y,t)\}`$ indeed satisfies (2), one needs to differentiate (25) with respect to $`x,y`$ and $`t`$ .

To be able to use this method we need to construct functions $`\mathrm{\Phi }_m`$, i.e. we need to solve for $`\mathrm{\Phi }_m`$ the system

$$\frac{_m\mathrm{\Phi }_n}{\mathrm{\Phi }_m\mathrm{\Phi }_n}=\frac{_mV_n}{V_mV_n},mn.$$
(27)
Observe that for $`mn`$

$$\frac{_nV_m}{V_mV_n}=\frac{\pi ^2\alpha _n}{\mathrm{sin}^2\pi (\gamma _m\gamma _n)};$$
(28)
this is a simple corollary of definitions of $`V_m`$ and $`W_m`$ and equations (20), (7).

Then the following functions satisfy equations (27):

$$\mathrm{\Phi }_m=\pi ^2_l\frac{H(\lambda )d\gamma }{\mathrm{sin}^2\pi (\gamma \gamma _m)},$$
(29)
where $`l`$ is an arbitrary closed contour on the branched covering $``$ such that its projection on the $`\lambda `$-plane does not depend on the branch points $`\{\lambda _n\}`$ and such that $`P_ml`$ for all $`m`$ ; $`H(\lambda )`$ is an arbitrary function on $`l`$ independent of $`\{\lambda _n\}`$. The proof of this fact is a simple calculation using equations (18), (12) and the link (13) between $`\gamma `$ and $`\mu .`$

Note that in this framework we can fix positions of branch points $`\{\lambda _{L+1},\mathrm{},\lambda _{2N2}\}`$ and consider the dependence of all functions on the remaining set of variables $`\{\lambda _1,\mathrm{},\lambda _L\},`$ where $`L2N2.`$

The following theorem summarizes our construction of solutions of the Boyer-Finley equation.

###### Theorem 1

Let functions $`\alpha _m(\{\lambda _n\})`$ and $`\gamma _m(\{\lambda _n\}),`$ $`m,n=1,\mathrm{},L`$, $`L2N2`$, be associated with an $`N`$-sheeted branched covering as described above. Define the potentials $`U(\{\lambda _n\})`$ and $`\psi (\{\lambda _n\})`$ to be solutions of the following system of equations:

$$\frac{U}{\lambda _m}=4\pi ^2\alpha _m,m=1,\mathrm{},L;$$
(30)
$$\frac{\psi }{\lambda _m}=\pi \alpha _m\mathrm{cot}\pi \gamma _m,m=1,\mathrm{},L.$$
(31)
Let the $`(x,y,t)`$dependence of branch points $`\lambda _n,n=1,\mathrm{},L`$, be governed by the following system of $`L`$ equations,

$$\pi ^2_l\frac{H(\lambda )d\gamma }{\mathrm{sin}^2\pi (\gamma \gamma _m)}=t+xV_m+yW_m,m=1,\mathrm{},L,$$
(32)
where

$$V_m=e^{2\pi i(\gamma _m+\psi )+U/2},W_m=e^{2\pi i(\gamma _m+\psi )+U/2};$$
(33)
$`l`$ is an arbitrary $`\{\lambda _n\}`$-independent contour on $``$ such that all $`P_ml`$; $`H(\lambda )`$ is an arbitrary summable $`\{\lambda _n\}`$-independent function on $`l`$.

Then the function $`U(\{\lambda _n(x,y,t)\})`$ satisfies the Boyer-Finley equation (1).

###### Remark 1

If an $`N`$-sheeted rational branched covering $``$ with two marked points $`Q_1,Q_2`$ is fixed, the solution of the Boyer-Finley equation constructed according to this theorem is defined by

(a) a functional parameter $`H(\lambda )`$ and

(b) a number $`L2N2`$, which has a meaning of the number of components $`\lambda _m`$ satisfying systems of hydrodynamic type (2) with characteristic speeds (33).

The application of theorem 1 in practice requires calculation of quadratures (30) and (31); besides that, one needs to resolve implicit relations (32) to find the dependence of $`\lambda _m`$ on $`(x,y,t)`$.

So far we were dealing with complex solutions of the Boyer-Finley equation (1); it is easy to formulate conditions on the parameters of our solutions which provide the reality of the function $`U`$.

Let us assume the function $`R(\gamma )`$ to satisfy the “reality condition”

$$\overline{R(\overline{\gamma })}=R(\gamma ).$$
(34)
Then the branch covering $``$ is invariant with respect to the antiholomorphic involution $`\tau `$, which acts on the points $`(\lambda ,\mu )`$ of the covering $``$ as follows:

$$\tau :(\lambda ,\mu )(\overline{\lambda },\overline{\mu }).$$
(35)
Assume also that both points $`Q_1,Q_2`$ are invariant with respect to $`\tau `$, i.e., $`\kappa _{1,2}`$. Let also the contour $`l`$ be invariant with respect to the involution and the function $`H(P)`$ satisfy the relation $`H(P)=\overline{H(P^\tau )}`$. Then one can choose the constants of integration in (30) and (31) such that the solution $`U(x,y,t)`$ of the Boyer-Finley equation given by theorem 1 is real.

Indeed, the invariance of the covering $``$ with respect to $`\tau `$ means that all $`\lambda _m`$ are either real or form conjugate pairs; the same holds for the set $`\{\mu _m\}`$. The expression (13) for the map $`\gamma `$ implies

$$\overline{\gamma (P^\tau )}=\gamma (P),$$
(36)
therefore, all $`\gamma _m\gamma (P_m)`$ are either imaginary, $`\overline{\gamma _m}=\gamma _m`$, or form anti-conjugate pairs. Applying complex conjugation to both sides of equation (18), we find that $`\alpha _m`$ are either real ($`\alpha _m`$ if $`\lambda _m`$) or form conjugate pairs, $`\alpha _m=\overline{\alpha _n}`$, if $`\lambda _m=\overline{\lambda _n}`$. This readily implies that one can choose the integration constant in the definition (30) of potential $`U`$ in such a way that $`U`$ is a real function of $`\{\lambda _m\}`$.

For completeness, we should also check that the reality condition does not contradict the solvability of system (32). Assume, for simplicity, that all $`\lambda _m`$ are real, i.e., all $`\gamma _m`$ are imaginary and all $`\alpha _m`$ are real. Then the potential function $`\psi `$ solving system (31) can be chosen to be imaginary, and both $`V_m`$ and $`W_m`$ (33) are real. Together with $`\overline{H(P^\tau )}=H(P)`$ and (36), it implies that both sides of equations (32) are real. Therefore, system (32) gives $`L`$ real equations for $`L`$ real variables $`\{\lambda _m\}`$, and generically has solutions.

Similar consideration applies when some $`\lambda _m`$’s form conjugate pairs; in this case the corresponding equations (32) will be conjugate to each other, and the number of real equations will again coincide with the number of real variables.

Acknowledgements We thank Tamara Grava and Yavuz Nutku for useful discussions.