text:
# Simulations of Galaxy Formation in a ΛCDM Universe II: The Fine Structure of Simulated Galactic Disks

## 1 Introduction

Some of the observed properties of galactic disks seem at odds with many of the “natural” trends expected in hierarchically clustering models and present a significant challenge to the current paradigm of structure formation on small scales (see, e.g., Sellwood & Kosowsky 2001, Navarro & Steinmetz 2000, Moore 2001). Qualitatively, the main difficulty lies in reconciling the early collapse and eventful merging history characteristic of the buildup of galactic dark matter halos with the many dynamical clues which point to a smooth assembly of the luminous component of galactic disks. This difficulty has made it necessary to postulate a substantial role for complex astrophysical processes in order to overcome some of these trends and to bring models of hierarchical assembly into agreement with observation.

The prime concern regards the fragility of centrifugally supported stellar disks to rapid fluctuations in the gravitational potential such as those stirred by mergers and satellite accretion events (Toth & Ostriker 1992, Quinn, Hernquist & Fullagar 1993, Velazquez & White 1999). As first computed in detail by Tóth & Ostriker (1992), the constraints on merger events undergone by a thin disk such as that of the Milky Way are very strict indeed. These authors find that less than $`10\%`$ of the disk mass within the solar circle could have been accreted in the past $`5`$ Gyr in the form of clumps, limiting severely the role of merging in the recent mass accretion history of the Milky Way.

Although these numbers might be relaxed somewhat by taking into account the coherent response of a self-gravitating disk to the accretion event (Huang & Carlberg 1997, Walker, Mihos & Hernquist 1996, Velázquez & White 1999), there is broad consensus that the presence of a dominant, thin, cold stellar disk points convincingly to a history of mass accretion where major mergers have played a minor role. Indeed, since the fragile circular orbits of disk stars are easily perturbed during mergers, the age of the oldest disk stars is often used to estimate the epoch of the last major merger. In the case of the Milky Way, stars as old as $`14`$ Gyr are found in the solar neighborhood, suggesting a protracted (and perhaps episodic) history of star formation in the disk (Wyse 2000, Rocha-Pinto et al 2000, Liu & Chaboyer 2000) and, consequently, a paucity of merger events at odds with typical merging histories of dark halos in cosmogonies such as the $`\mathrm{\Lambda }`$CDM model.

The need for ordered and smooth settling of gas into thin galactic disks leads to a few basic predictions of hierarchical scenarios that can be contrasted directly with observation. For example, because angular momentum results from torques operating before collapse, it correlates directly with the radius and time of turnaround (Navarro & Steinmetz 1997). Thus, material that turns around later is expected to have higher angular momentum, with two major implications: (i) disks are expected to be physically smaller in the past, compared with systems of similar rotation speed identified at present, and (ii) disks form from the inside out, as higher angular momentum material, which collapses later, should settle preferentially in the outer regions of the disk (Mo, Mao & White 1998).

Although overall these expectations are not grossly inconsistent with current data, a number of well established observations are worryingly difficult to accommodate within this general scenario. One of these worries concerns the origin of the complex vertical structure of galactic disks. Current datasets suggest that most (if not all) galactic disks are built out of two major dynamical components, referred to generally as the “thick” and “thin” disks (Dalcanton & Bernstein 2002 and references therein). Although varying in importance from galaxy to galaxy, the thick disk appears to be old ($`>10`$-$`12`$ Gyr in the case of the Milky Way, Gilmore, Wyse & Jones 1995), poor in metals (Prochaska et al 2000), and of similar radial extent as the thin disk component (Wyse 2000).

These trends have led to the popular conjecture that the Milky Way’s thick disk has its origin in an early thin disk of velocity and size comparable to today’s but “thickened” by the accretion of a satellite $`10`$-$`12`$ Gyr ago (Gilmore, Wyse & Norris 2002). Appealing as this idea may be from a dynamical standpoint, it is rather unattractive in a hierarchical scenario, where the mere existence of such large, fast-rotating disks at $`z2`$ (corresponding to roughly $`12`$ Gyr ago in the $`\mathrm{\Lambda }`$CDM scenario) counters the “natural” prediction of the model. A related difficulty concerns the presence of old disk stars in the vicinity of the Sun, which might be taken to imply that even at high redshift the thin disk already extended as far as the solar circle and has remained essentially undisturbed dynamically since. Such a long period of quiet dynamical evolution is difficult to reconcile with the active merging expected at early times in the $`\mathrm{\Lambda }`$CDM cosmogony, and perhaps even with the presence of the Galaxy’s spheroidal component, which is typically ascribed to a relatively recent major merger event.

The discussion above illustrates the still unsettled account of the origin of the Galactic disk(s) in the hierarchical formation paradigm, and draws notice to the fundamental role that accretion events, be they mergers or satellite disruptions, have played in determining the structure of galactic disks. In this paper—the second in a series analyzing the dynamical properties of galaxies simulated in the “concordance” $`\mathrm{\Lambda }`$CDM paradigm—we examine these issues using a numerical simulation with high numerical resolution. The simulation is the same as described by Abadi et al (2002) in Paper I of this series; that paper confronts the global dynamical and photometric properties of the simulated galaxy with normal spirals, while the present one concentrates on the multi-component nature of the stellar disk and on its dynamical origin.

The paper is organized as follows. In section 2 we briefly discuss, for completeness, numerical details of the simulation; § 3 discusses the identification of distinct dynamical components in the simulated galaxy, as well as the implications of these results in the context of the formation of the thick and thin disks of the Milky Way. Finally, § 4 summarizes our main conclusions.

## 2 The Simulation

A complete description of the numerical simulation analyzed in this paper is presented in Paper I (Abadi et al 2002), where the interested reader may find full details. For completeness, we provide here a brief description of the numerical method and of the simulation.

### 2.1 The Code

The simulation described here was performed with GRAPESPH, a particle-based, fully three-dimensional Lagrangian hydrodynamical code that combines the flexibility and adaptability of the Smoothed Particle Hydrodynamics technique with the speed of the special-purpose hardware GRAPE for computing gravitational interactions between particles (Steinmetz 1996). The version used here includes the self-gravity of gas, stars, and dark matter, hydrodynamical pressure and shocks, Compton and radiative cooling, as well as the heating effects of a photoionizing UV background (see Navarro & Steinmetz 1997, Steinmetz & Navarro 1999, for more details).

Star formation is handled in GRAPESPH by means of a simple recipe for transforming gas particles in collapsing, dense regions into stars at a rate roughly proportional to the local dynamical timescale of the gas. After formation, star particles are only affected by gravitational forces, but they devolve energy and mass to their surroundings in a crude attempt to mimic the energetic feedback from evolving stars and supernovae. A fraction, $`ϵ_v=0.05`$, of this energy is invested in raising the kinetic energy of gas neighboring active star forming regions. These motions are dissipated by shocks on a timescale longer than the local cooling time, allowing for lower star formation efficiencies and longer effective timescales for the conversion of gas into stars. Still, our choice of parameters results in just a minor fraction of cooled, dense gas being reheated and returned to intergalactic space in diffuse form. The star formation history in the simulated galaxy thus roughly traces the rate at which gas cools and collapses to the center of dark matter halos.

The mass devolved by star particles during feedback events is added to gas particles surrounding young stars, and tracked in the code as a crude measure of the metal content of the gas. This “metallicity” is in turn inherited by the stars each gas particle spawns. By all accounts, this treatment of metal enrichment is quite rudimentary and, therefore, we restrict the use of stellar metallicities to the modeling of the spectrophotometric contribution of each star particle to the luminous output of a galaxy, as described by Contardo, Steinmetz & Fritze-von Alvensleben (1998).

### 2.2 The Simulation

The simulation evolves a region that forms, at $`z=0`$, a galaxy-sized dark matter halo in the low-density, flat Cold Dark Matter ($`\mathrm{\Lambda }`$CDM) scenario (Bahcall, Ostriker & Steinhardt 1999). This model is fully specified by the following choice of cosmological parameters111 We express the present value of Hubble’s constant as $`H(z=0)=H_0=100h`$ km s-1 Mpc-1: $`\mathrm{\Omega }_0=0.3`$, $`h=0.65`$, $`\mathrm{\Omega }_\mathrm{b}=0.019h^2`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, where $`\mathrm{\Omega }_0`$, $`\mathrm{\Omega }_b`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ represent the present-day contribution to the mass-energy density of matter, baryons, and a cosmological constant, respectively, in units of the closure density. The power spectrum is normalized so that at $`z=0`$ the linear rms amplitude of mass fluctuations in $`8h^1`$ Mpc spheres is $`\sigma _8=0.9`$, and we assume that there is no “tilt” in the initial CDM power spectrum.

At $`z=0`$ the dark matter halo under consideration has a circular velocity of $`V_{200}=134`$ km/s and a total mass of $`M_{200}=5.6\times 10^{11}h^1M_{}`$, measured at the virial radius, $`r_{200}=134h^1`$ kpc, where the mean inner density contrast (relative to the critical density for closure) is $`200`$. This region was identified in a cosmological simulation of a large periodic box ($`32.5h^1`$ Mpc on a side) and resimulated at higher resolution, including the tidal field of the original simulation. The high-resolution region (an “amoeba”-shaped region contained within a cube of $`3.4h^1`$ comoving Mpc on a side) is filled at the initial redshift, $`z_i=50`$, with the same number of gas and dark matter particles. The gas and dark matter particle mass is $`m_g=2.1\times 10^6h^1M_{}`$ and $`m_{dm}=1.2\times 10^7h^1M_{}`$, respectively. We adopt a Plummer softening scalelength of $`0.5`$ (physical) kpc for all gravitational interactions between pairs of particles.

The baryonic mass of the final galaxy is $`10^{11}M_{}`$, equivalent to roughly $`36,000`$ gas particles (or more than $`100,000`$ star particles). This resolution enables a detailed study of the dynamical and photometric properties of the simulated galaxy, including the identification of different populations of stars according to age and/or kinematics. Galaxy luminosities are computed by simply adding up the luminosities of each star particle, taking into account the time of creation of each particle (its “age”) and its metallicity, as described in detail by Contardo, Steinmetz & Fritze-von Alvensleben (1998). Corrections due to internal absorption and inclination are neglected, except for a temporary wavelength-dependent dimming introduced to take into account the gradual dispersal of the dust clouds that obscure the formation sites of young stars (Charlot & Fall 2000). Further details may be consulted in Paper I (Abadi et al 2002).

## 3 Results

The large-scale formation of the simulated galaxy is shown in Figures 1-3 of Paper I, and follows the typical assembly process of dark matter halos in the $`\mathrm{\Lambda }`$CDM cosmogony. The region destined to form the halo collapses early ($`z<5`$) into a pancake-like structure criss-crossed by filaments along which matter drains into collapsed dark matter halos. Gas cools and condenses at the center of these halos, turning swiftly into stars as it settles at the bottom of these potential wells. These progenitors gradually merge into a single system at $`z=0`$ which is the subject of our analysis here.

The mass accretion history is punctuated by a succession of major mergers at high redshift, although by $`z1`$ the merging activity of the galaxy is over and the dark halo evolves little thereafter. The cooled (and cooling) gas within the halo, on the other hand, keeps flowing towards the center and settles gradually into a dynamically cold, centrifugally supported disk. The stellar component of the central galaxy increases by $`40\%`$ in mass since $`z1`$, and is shown at $`z=0`$ in the top-left panel of Figure 1. Much of the late accreting material is clearly visible in the final galaxy as a well-defined disk of young stars surrounding a more centrally concentrated spheroid of older stars. Photometrically, this structure is reminiscent of that of early-type spirals and was analyzed in detail in Paper I. We concentrate here on the detailed dynamical analysis of the galaxy at $`z=0`$ and on the origin of its major components.

### 3.1 Dynamical Decomposition

The existence of a well defined, thin disk in the stellar and gaseous components of the galaxy (see top-left panel of Figure 1) defines unambiguously a symmetry axis which we shall hereafter call the “$`z`$-axis”, for short. Figure 2 shows $`J_z`$, the $`z`$-component of the specific angular momentum of all stars within the “luminous” radius of the central galaxy, $`r_{\mathrm{lum}}=21`$ kpc, plotted versus the specific binding energy of each star, $`E`$. Binding energies are computed relative to the total mass of the system within the virial radius, $`r_{200}`$. Particles with positive or negative $`J_z`$ circulate around the center of the galaxy in different directions; we choose $`J_z`$ to be positive so that most stars are “co-rotating” in Figure 2.

The available parameter space in the $`J_z`$-$`E`$ plane is delineated by the maximum and minimum angular momentum consistent with a given value of the specific energy: those of the co- and counter-rotating circular orbits, respectively. The co-rotating circular orbits are shown by the solid curve labelled $`J_{\mathrm{circ}}(E)`$ in Figure 2. The gas particles within $`r_{\mathrm{lum}}`$ are on a thin, centrifugally supported disk, and follow closely the $`J_{\mathrm{circ}}(E)`$ curve, as shown by the magenta dots in Figure 2. The prevalence of circular orbits in the gaseous disk can also be seen in the distribution of the “orbital circularity” parameter, $`ϵ_J=J_z/J_{\mathrm{circ}}(E)`$, which is sharply peaked about unity, as shown in the top inset of Figure 2.

The circularity distribution of stars also shows a sharp peak at $`ϵ_J1`$, but is quite broad and has approximately $`25\%`$ of stars on counter-rotating orbits (see solid line labelled “all stars” in the bottom inset of Figure 2). The shape of the distribution suggests the presence of at least two distinct dynamical components; a hot spheroid with as many stars in co- as in counter-rotating orbits, and a cold disk of stars with circularities resembling those of the gas and sharply peaked about unity. The curves labelled “spheroid” (long-dashed) and “thin disk” (dotted) indicate a possible decomposition into such components. This decomposition assumes that the spheroid has little net rotation, and that the circularity distribution of the thin disk is similar to that of a cold disk such as the thin disk of the Milky Way. Indeed, the circularity distribution of the “thin disk” identified in the bottom inset of Figure 2 is similar to that of a $`\sigma /V_{\mathrm{rot}}0.2`$ disk, where $`\sigma `$ and $`V_{\mathrm{rot}}`$ are the velocity dispersion and the rotation speed of the disk, respectively.

The top and bottom right panels of Figure 1 show the spatial distribution of stars selected according to the circularity distributions labelled “spheroid” and “thin disk” in Figure 2, respectively. These panels confirm the intuitive expectation that the spheroid ought to be roughly spherical in shape whereas the thin disk is highly flattened, with an aspect ratio exceeding $`10`$:$`1`$.

The decomposition into spheroid and thin disk proposed here leaves, however, a number of stars on prograde orbits unassigned to either component (see the short-dashed curve labelled “thick disk” in Figure 2). These stars share the same sense of rotation as the thin disk but rely less on rotation for their support; as a result, they are distributed spatially on a “thick disk” shown in the bottom left panel of Figure 1.

It must be recognized that the identification of three distinct dynamical components, although compelling in light of the shape of the stellar circularity distribution, is nevertheless somewhat arbitrary. One may question, for example, our assumption that the spheroid does not rotate, or our assignment of essentially all stars on circular orbits to the thin disk (realistic models for the thick disk would typically have a sprinkling of such orbits as well). We have experimented with several possible decomposition procedures and report here results which, although they depend quantitatively on the exact definition adopted for each component, are qualitatively independent of such details and hold for every plausible decomposition procedure. We list the main parameters of each component in Table 1 and concentrate below on their origin as well as on their imprint on the photometric properties of the galaxy.

### 3.2 Photometric vs Dynamical Components

As discussed in Paper I, the surface brightness profile of the simulated galaxy resembles that of early-type spirals, and can be described very well by combining a central $`R^{1/4}`$ spheroid with a surrounding exponential disk. This is shown by the solid circles in Figure 3, where we show the $`I`$-band surface brightness profile of the simulated galaxy seen face-on. The solid line through the circles corresponds to the best $`R^{1/4}`$+exponential fit, whereas the curves labelled “$`R^{1/4}`$” and “exp” indicate the contribution of each of these components to the fit. The best $`R^{1/4}+`$exponential fit splits the total light in the $`I`$ band roughly equally between the disk and the spheroid, as in early type spirals such as UGC615. The virtues and shortcomings of this comparison are discussed in detail in Paper I. Here we focus on the dynamical interpretation of the fitting procedure.

Does the photometric decomposition describe accurately the true dynamical importance of each component? Is the disk truly exponential? Is the spheroid really well approximated by de Vaucouleurs’ law? Our simulation offers intriguing clues to these questions. The inverted triangles and filled squares in Figure 3 show the contribution of the dynamically-identified thin disk and spheroid to the face-on $`I`$-band surface brightness profile, respectively.

The spheroid contributes most of the $`I`$-band light ($`60\%`$); the thin disk $`30\%`$ and the thick disk (shown with upright triangles in Figure 3) the remaining $`10\%`$. The total contribution of the (thick+thin) disk thus is of order $`40\%`$, similar to that inferred from the photometric decomposition. This is somewhat surprising, since the individual profiles deviate strongly from the assumed $`R^{1/4}`$ and exponential laws. Indeed, neither component can be adequately approximated by one such law; both show a roughly exponential decay in the outer regions as well as a sharp upturn near the center well described by a de Vaucouleurs’ profile.

The contribution of the thin disk to the total stellar mass of the system is only $`17\%`$; its prominence in the surface brightness profile is due to its younger age (and thus lower mass-to-light ratio) compared with the spheroid. This is shown in Figure 4; the mean luminosity-(mass) weighted age of the thin disk is $`3.6`$ ($`5.3`$) Gyr; a full $`7.3`$ ($`6.0`$) Gyr younger than the spheroid. For a Scalo IMF (used throughout this analysis), this age difference results in the $`I`$-band mass-to-light ratio of the thin disk being $`2.2`$ times brighter than the spheroid, accounting for the importance of the disk in the photometric decomposition. The thick disk is of intermediate age between the thin disk and the spheroid; its luminosity (mass)-weighted age is $`9.8`$ ($`10.5`$) Gyr.

Clearly, inferring the structural properties of the disk and spheroidal components from photometric decomposition alone is not straightforward. This may have important implications for studies that attempt to assess the relative importance of cold, centrifugally-supported disks and of hot, “pressure”-supported spheroids through photometric decomposition of large samples of galaxies (Schechter & Dressler 1987, Benson, Frenk & Sharples, 2002). In particular, our study shows that it is critical to account carefully for the age difference between components and its effect on the mass-to-light ratios: assuming similar mass-to-light ratios for the disk and spheroid would result in dramatic overestimation of the dynamical importance of the disk component.

The results shown in Figure 3 have repercussions as well for studies that rely on photometric decomposition to gain insight into the process of bulge formation. For example, the tight link between the spatial scalelengths of the spheroidal and disk components derived photometrically might not necessarily reflect a physical association between the two, as is usually interpreted in models that view the formation of the bulge as a result of secular evolution in the disk (MacArthur, Courteau, & Holtzman 2002, Courteau, de Jong, & Broeils 1996). For example, the half-mass radii of the thin disk and spheroid are $`5.0`$ and $`0.7`$ kpc, respectively, which should be contrasted with the $`7.8`$ and $`1.1`$ kpc that would be obtained from the best $`I`$-band $`R^{1/4}`$+exponential fit assuming constant mass-to-light ratios for each component.

Inferring the radial extent of the disk from photometric data alone thus leads to a $`50\%`$ overestimate of the half-mass radius. This is an important correction, since disk scalelengths estimated photometrically are usually used to estimate the angular momentum of galaxy disks (see Paper I and references therein). It is important that angular momenta computed from numerical simulations are estimated in the same way as observations in order to ensure that any discrepancies reflect genuine differences in angular momentum content and not systematic uncertainties in the estimation procedure.

With only one example, it would be premature to validate (or invalidate) the conclusions of the studies mentioned above, but clearly caution is highly recommended when using photometric data alone to assess dynamical theories for the origin of the components of the galaxy. We intend to return to this question when we have completed a more statistically significant sample of simulated galaxies.

### 3.3 Dynamical Components and Stellar Populations

The dynamical distinction at the basis of the decomposition procedure described in § 3.1 is closely linked to the time and mode of assembly of each component. As shown in Figure 4a, the spheroid is old and very well mixed dynamically; its mass-weighted and luminosity-weighted ages are quite uniform across the galaxy, with little sign of a radial gradient. The thin disk, on the other hand, is fairly young and has a very well defined radial gradient as a result of its inside-out assembly; the luminosity (mass) weighted age drops from $`5`$ ($`6`$) Gyr at $`R1`$ kpc to just about $`2`$ ($`3`$) Gyr at $`R=20`$ kpc. The age difference between components is reflected in their colors, as shown in Figure 4b. The color profile of the simulated galaxy as a whole (solid circles) shows a mild gradient, becoming bluer farther from the center, in good agreement with that of the Sab galaxy UGC615 (open circles). This trend is due almost entirely to the age gradient in the thin disk component shown in panel (a) of Figure 4; metallicity plays only a minor role in this result.

Our results confirm two of the basic expectations of galaxy models based on hierarchical accretion. (i) Spheroids are predominantly old and very well mixed dynamically by merger events which weaken any preexisting age or metallicity gradient in the progenitors. This is, of course, consistent with the Searle & Zinn (1978) interpretation of the lack of obvious gradients in the stellar halo of the Milky Way as due primarily to its formation through merging of protogalactic fragments. (ii) Centrifugally-supported disks, on the other hand, are built from the inside out and are on average much younger than spheroids. These results seem broadly consistent with observational data for external galaxies and for our own Milky Way. We turn next to the detailed dynamical origin of each component and to the subtle but important distinction between the thin and the thick disk components.

### 3.4 The Thin and Thick Disks

We have identified a “thick disk” component with stars on co-rotating (but rather eccentric) orbits which do not belong to either the thin disk or the spheroid, according to the decomposition procedure outlined in § 3.1. The combined thin+thick disk component (hereafter “the disk”, for short) resembles that of the Milky Way in several respects. Consider, for example, the vertical distribution of disk stars, shown in Figure 5. The $`z`$-density profile deviates clearly from a simple exponential law, but is well approximated by the combination of two exponentials, with $`575`$ pc and $`2.7`$ kpc scaleheights respectively (see filled circles and solid lines in Figure 5).

Interestingly, stars identified dynamically as belonging to the thin disk follow closely an exponential profile as well, with scaleheight $`570`$ pc. The density profile of the thick disk is not exactly exponential, but the median height of its stars is not very different from that inferred from the double-exponential fit. The disk rotates significantly slower at higher $`z`$, as shown in the bottom panel of Figure 5. The mean rotation speed in the outskirts of the disk ($`10`$ kpc ¡ R ¡ $`15`$ kpc, where the rotation curve is approximately flat, see Paper I) drops from $`240`$ km/s on the plane of the disk to $`160`$ km/s at $`z3`$ kpc.

Are the thick and thin disks really distinct components, or is the double exponential profile shown in Figure 5 just a manifestation of a rotationally-supported structure intrinsically more complex than a simple exponential? To be meaningful, a dynamical decomposition between thick and thin disk components must also identify other independent attributes, such as age or metallicity, which neatly divide the stellar population into similar components. In the case of the Milky Way, the metallicity and age distributions of disk stars support the idea that the thin and thick disks are truly distinct components of independent origin. Besides having higher velocity dispersion, thick disk stars are typically older and more metal poor than thin disk stars, and appear to divide “naturally” into two components in the age/metallicity plane (see, e.g., §10.4.3 of Binney & Merrifield 1998).

The distinction between thin and thick disks is also often illustrated by simultaneously ranking stars by their metal abundance and by their vertical speed (see Freeman 1991). Since our treatment of metal enrichment is rather rudimentary (§ 2.1), we illustrate the same idea using age rather than metallicity as one of the ranking variables. The break at $`8`$ Gyr seen in the top panel of Figure 6 is highly suggestive of the presence of two components of distinct nature and origin: one young and cold, well approximated by a velocity dispersion of $`\sigma _z30`$ km s-1, overlapping an older population of stars with a velocity dispersion roughly three times as high. (A Gaussian distribution of velocities would show as a straight line in this plot.)

This interpretation is confirmed by the velocity dispersion-age relation for disk stars shown in the bottom panel of Figure 6, where a sharp upturn in velocity dispersion is noticeable for stars older than $`8`$ Gyr. These results support the view that the thin and thick disks are truly distinct components, and point to a major occurrence in the assembly of the galaxy roughly $`8`$ Gyr ago: this timescale coincides with the epoch of the last major accretion event.

Major accretion events usually trigger profuse star formation and thus show clearly in the distribution of stellar ages shown in Figure 7. For example, the sharp peak at age $`\tau 9.5`$ Gyr corresponds to a merger event which increases the stellar mass of the galaxy by almost $`50\%`$. Similarly, the broad peak at $`\tau 8`$ Gyr is associated with the nearly simultaneous accretion of two satellites; these and the associated starburst boost the stellar mass of the galaxy by $`20\%`$. Subsequently, the main galaxy is disturbed only by the accretion of two minor satellites (at $`z0.7`$), each contributing less than $`6\%`$ of the current stellar mass of the galaxy. At lookback times less than $`8`$ Gyr most stars form (and remain) in a thin disk as a result of the gradual transformation into stars of cooled gas accreted smoothly into the main body of the galaxy.

Clearly, the sharp decline in merger activity which occurred $`8`$ Gyr ago is the watershed event responsible for initiating the formation of the thin disk and for establishing the neat separation between the thin and thick disk components. The bulk of thick disk stars is therefore older than $`8`$ Gyr and the majority of thin disk stars form after that epoch. We explore next the origin of the thick disk in the simulated galaxy and the applicability of these findings to the formation of the Milky Way’s thick disk.

### 3.5 The Origin of the Thick Disk

As discussed in § 1, dynamical data in the solar neighborhood are usually interpreted within a scenario where the thick disk results from the tidal heating of an early thin stellar disk. Does the thick disk in the simulated galaxy form in this way? We can test this proposal directly by tracing the stars assigned to the thick disk at $`z=0`$ and by inspecting their configuration at the time of the last major merger.

The result of this exercise is summarized in the bottom panel of Figure 8. This figure shows, as a function of age, the mass fraction of each component which formed within the main progenitor of the simulated galaxy. The evolution of this “in situ”-formed mass fraction is similar for the thin and thick disks, and indicate that essentially all disk stars younger than $`8`$ Gyr were formed within the main galaxy, a result consistent with the lack of major accretion events since $`z1`$.

On the other hand, the majority of old disk stars were formed in systems other than the main progenitor; only about $`10\%`$ of disk stars older than $`10`$ Gyr were formed in the main progenitor of the galaxy. This implies that the bulk of old disk stars are actually tidal debris from disrupted satellites. Since the thick disk is quite old (only $`12\%`$ of its stars are younger than $`8`$ Gyr, see the top panel of Figure 8), this rules out the hypothesis that the thick disk is the remnant of a stirred pre-existing thin disk and points to tidal disruption as the main process of assembly of the thick disk. About $`70\%`$ of thick disk stars older than $`8`$ Gyr are debris from tidally disrupted satellites; in total, more than $`60\%`$ of thick disk stars of all ages were brought in by satellites, the rest were mostly formed during the accretion events themselves.

Interestingly, although the bulk of the thick disk consists of debris from tidally disrupted satellites, the fraction of spheroid stars identified as tidal debris is only $`50\%`$ of the total. The other half of the stars in the spheroid formed in a disk-like component in the main galaxy but were dispersed into the spheroid as a result of a major ($`1`$:$`2`$) merger at $`z1.5`$ ($`9.5`$ Gyr ago). Thus, our results suggest that, contrary to what might be naively expected, it would be easier to find remnants of tidal stripping events in the old disk component rather than in the spheroid.

Considering the thick disk as the accumulated tidal debris of disrupted satellites rather than as an early thin disk thickened by an accretion event helps to resolve the angular momentum puzzle alluded to in § 1. The high angular momentum content of the thick disk does not reflect the presence of a large disk at high redshift but rather the orbital angular momentum of disrupted satellites, some of which were aggregated into the main galaxy relatively recently. In our simulation, for example, half of the stars in the thick disk are older than $`10`$ Gyr, but were fully accreted into the main progenitor only $`6`$ Gyr ago. We explore below how this result also affects the interpretation of very old stars on circular orbits in the outskirts of the disk galaxies.

### 3.6 The Old Thin Disk

As discussed in § 1, the existence of disk stars which predate the epoch of the last major merger but which are nevertheless found today on nearly circular orbits is puzzling in hierarchical assembly models of disk formation. Debris from tidally disrupted satellites offers again a plausible explanation for the presence of such a component. The simulated galaxy has a well defined “old thin disk” population; indeed, Figure 8 shows that $`15\%`$ of the thin disk stars are older than $`10`$ Gyr, and therefore formed long before the time of the last major accretion event. As shown in the bottom panel of the same Figure, essentially all of these stars were formed in satellites and brought into the main galaxy much later; the last of these “old thin disk” stars was only incorporated into the main galaxy roughly $`6`$ Gyr ago.

Figure 8 thus demonstrates that tidal debris from satellites can contribute not only to the spheroid (as is usually assumed) but also to the thin and thick disk components. A satellite’s contribution to either component (or to the spheroid) depends mainly on its orbit and on the degree to which dynamical friction circularizes the orbit before disruption. In order to contribute significantly to the disk, a satellite must have an orbital plane roughly coincident with that of the disk and be dense enough to survive disruption until its orbit has been circularized within the disk.

This intuitive expectation is confirmed by the data presented in Table 2, where we list the structural properties and orbital parameters of four different satellites accreted after $`z1`$. This Table lists the satellite’s orbital circularity, as well as the fraction of stars that end up, at $`z=0`$, in the spheroid, thin disk, and thick disk components, respectively. Since the orbital properties change significantly with time, we have chosen to measure them at $`z_{\mathrm{acc}}`$, the time of their last apocentric passage before significant disruption takes place.

As expected, satellites on orbits of high circularity are seen to contribute principally to the disk component rather than to the spheroid, and viceversa. Closer inspection shows that it is the satellite’s core that contributes primarily to the thin disk whilst stars in its outskirts end up primarily in the thick disk or spheroid (see Figure 9). This is fully consistent with the notion that dynamical friction circularizes the orbit of the satellite prior to disruption: stars at the core of the satellite are more resilient to disruption and therefore end up on more circular orbits than those stripped earlier during the satellite decay process.

The discussion above implies that, although the disk component is formed overall from the inside out (see, e.g., Figure 4), a small fraction of old stars on circular orbits is expected at essentially all radii as a result of satellite accretion. Is this a result of general applicability or just a fluke in this simulation? The answer depends on the incidence of capture of satellites on orbital planes roughly coincident with that of the disk of the galaxy.

With the obvious caveat that it is perilous to generalize from a single example, there is indication that the accretion of such satellites might not be uncommon in hierarchical scenarios. As Table 2 indicates, all $`4`$ satellites accreted by the galaxy since $`z1`$ co-rotate with the disk and, after disruption, contribute roughly $`40\%`$ of their stars to the disk. This result, as well as the relatively high frequency of satellites with orbital planes coincident with the disk, are apparently related to the fact that much of the late accreting material (which effectively forms the disk) flows in through a filament.

This filament inevitably contains previously collapsed, condensed subsystems; satellites whose orbits share the properties of the material responsible for the formation of the disk. If dense enough, these satellites survive disruption until after their orbits are circularized and contribute a substantial fraction of their stars to the disk. We expect this to be a common feature of the hierarchical formation process of a disk galaxy, implying that the majority of galactic disks ought to have a significant population of old stars on orbits similar to those of young disk stars but with distinct properties such as age and metallicity. We intend to explore fully the consequences of this idea in future papers of this series. We discuss briefly below the implications of these findings for the Milky Way galaxy.

### 3.7 Implications for the Milky Way

The discussion above illustrates the important clues to the accretion history of the galaxy retained in the detailed dynamics and structure of the oldest stellar populations. In particular, it shows conclusively that relicts of past accretion and disruption events might be found today not only in the stellar halo of the Galaxy, as is commonly assumed (Helmi et al. 1999, Majewski et al 2000), but also in the cold, centrifugally supported disk component. As mentioned above, tidal debris is proportionally more abundant in the old disk component of the simulated galaxy than in the old spheroid; $`85\%`$ of disk stars older than $`10`$ Gyr are the relicts of disrupted satellites, well in excess of the $`55\%`$ of the old stellar spheroid contributed by satellite debris. This suggests that it would be profitable to search for signatures of accretion events in the disk as well as in the halo components; indeed, such relicts may have already been detected in the unexpected dynamics of stars above and below the plane of the Galaxy reported by Gilmore et al (2002).

Could the “old disk” in the Milky Way be made up of the remains of disrupted satellites? It is important to verify that this proposal is consistent with the observation that the radial distribution of old disk stars is not very different from that of the younger thin disk component. This is especially true because stars from a single disrupted satellite are not uniformly mixed radially but rather form a torus-like structure whose radius is determined by the location at which final disruption occurs. This is clearly shown in Figure 9, where the stars of satellite $`1`$ (see Table 2) are shown at different times during its disruption. Once the disruption is complete, the satellite stars are distributed along a thin “ring” of stars at a radius $`1.7`$ kpc from the center.

On the other hand, several satellites contribute to the old disk component, and each contributes in varying amounts at different radii. This is shown in the left panel of Figure 11; the debris from satellite $`1`$ dominates the inner regions ($`R<5`$ kpc) of the old thick disk, whereas debris from satellite $`2`$ is the principal contributor to the outer regions. Satellites $`3`$ and $`4`$ contribute fewer stars at all radii because their orbits are far less circular (see Table 2).

Fewer satellites are likely to contribute to the thin disk than to the thick disk. This implies that, at given radius, the old thin disk might be dominated by the contribution of a single satellite. This is indeed the case in the simulated galaxy at $`R<10`$ kpc, as shown in Figure 11b. Systematic variations in the properties of old disk stars as a function of their rotational support are thus expected and should contain intriguing clues to the early assembly process of the disk. This subtle modulation might be observable in our own Milky Way, once surveys such as those planned for the space-based GAIA mission (Perryman et al 2001) and the ground-based RAVE project (Steinmetz 2002)—where full dynamical information will be available for large samples of stars well beyond the solar neighborhood—are completed.

Overall, the radial overlap of debris from several disruption events results in a radial distribution of old disk particles that differs little from that of the young disk component, in agreement with observation. This is illustrated in Figure 12, where we show the surface density profile of the “young” and “old” disk components. There are no major differences in the radial distribution of the old and young disks; the half-mass radius of the old thin (thick) disk is $`3.9`$ ($`5.6`$) kpc, comparable to that of the young thin (thick) disk, $`5.2`$ ($`3.9`$) kpc, despite the large difference in the mean age of each component. The late accretion of satellites with orbits coincident with the disk thus explains the puzzle alluded to in § 1 regarding the radial extent (and angular momentum content) of the older, thicker disk component of the Galaxy.

If much of the old thick (and thin) disk is indeed debris from tidally disrupted satellites, their metallicities must reflect the metallicity of the satellites themselves. Since metallicity and luminosity are strongly correlated in dwarf systems, then one may use the metallicity of the old disk component to infer the luminosity of the disrupted satellite. For example, if the metallicity of the thick disk is taken to be \[Fe/H\]$`0.6`$, one would infer that the luminosity of its progenitor satellite had a luminosity comparable to that of the LMC, $`M_V18.5`$ (see Figure 8 of van den Bergh 1999).

Another important implication of this idea is that the same satellite can contribute both to the old thin and thick disks. This implies that these two components should exhibit a clear overlap in properties, reflecting their common origin in a single system. This is reminiscent of the results of Wyse & Gilmore (1995), who argue that the dynamical and metallicity distributions of disk stars can only be explained if stars with abundances \[Fe/H\] $`<0.4`$ contribute in substantial amounts to both the thick and thin disk components.

We end this discussion with the obvious caveat that, although we find the formation mechanism of the old disk proposed here attractive and applicable to the Milky Way, the simulated galaxy is not intended to be a model of the Galaxy. Further simulation work is required to ascertain the general validity of the results reported here and their applicability to the Milky Way and other spirals.

## 4 Summary

We present a detailed analysis of the dynamical components of a galaxy simulated in the “concordance” $`\mathrm{\Lambda }`$CDM cosmogony. The galaxy forms in a dark matter halo chosen so that mergers and accretion events are unimportant dynamically after $`z1`$. Star formation and feedback parameters are such that the star formation history of the galaxy is largely driven by the rate at which gas cools and collapses within dark halos: this is a conservative choice where feedback effects are relatively unimportant. The shortcomings of this assumption concerning the global photometric and structural properties of the simulated galaxy are discussed in Paper I. We focus here on the multi-component nature of the stellar disk and on their dynamical origin. Our main results can be summarized as follows.

* At $`z=0`$, two separate stellar components are easily distinguishable solely on the basis of the orbital parameters of stars in the galaxy: a slowly rotating, centrally concentrated spheroid and an extended disk-like component largely supported by rotation.
* These components are also recognized in the surface brightness profile of the galaxy, which can be very well approximated by the superposition of an $`R^{1/4}`$ spheroid and an exponential disk, in agreement with observations of early-type spirals.
* Neither does the dynamically identified spheroid follow closely de Vaucouleurs’ law nor is the disk purely exponential, a result which calls for caution in the interpretation of estimates of the dynamical importance of the disk and spheroid in traditional photometric decomposition techniques.
* The spheroid is old, and has essentially no stars younger than the last major accretion episode; $`8`$ Gyr ago for the system we consider here. The majority of thin disk stars, on the other hand, form after the merging activity is over and have a mean age of $`5`$ Gyr.
* The disk may be further decomposed into two well defined subcomponents: a thin, dynamically cold disk of stars on nearly circular orbits and a thicker disk with orbital parameters transitional between the thin disk and the spheroid. Supporting evidence for the true presence of these two distinct components is found, as in the Milky Way, in the double-exponential vertical structure of the disk and in abrupt changes in the vertical velocity distribution as a function of age.
* The bulk ($`60\%`$) of the thick disk consists of the tidal debris of satellites whose orbital plane was coincident with the disk and whose orbits were circularized by dynamical friction prior to full disruption. This trend becomes more pronounced with age; $`90\%`$ of stars older than $`10`$ Gyr and presently in the thick disk component were brought into the disk by satellites.
* A significant fraction ($`15\%`$) of thin-disk stars are old enough ($`>\mathrm{\hspace{0.17em}10}`$ Gyr) to predate the last major accretion event. The bulk of this unexpected population of old stars on nearly circular orbits share a common origin with the old thick disk: they are the remains of the cores of disrupted satellites. Interestingly, only one in two stars belonging to the old spheroid are tidal debris; the rest may be traced to a major merger event that disperses the luminous progenitor at $`z1.5`$ and seeds the formation of the spheroid.

Our results offer clues to understanding a number of observational trends that challenge the standard hierarchical disk assembly process. In particular, the presence of an old disk component made up primarily of tidal debris from disrupted satellites helps to explain: (i) the presence of a significant number of old stars on circular orbits in the outskirts of galaxies like the Milky Way, and (ii) why the specific angular momentum and radial extent of the thick and thin disk components are comparable in spite of the significantly different ages of their individual stars. It also offers a “natural” explanation for the clear dynamical and evolutionary distinction between the thin and thick disk components: the thick disk is mostly tidal debris from disrupted satellites whereas the young thin disk consists mostly of stars formed “in situ” after the merging activity abates.

These findings highlight the role of satellite accretion events in shaping the disk—as well as the spheroidal—components of a galaxy and reveal some of the clues to the assembly process of the galaxy preserved within the detailed dynamics of old stellar populations. We emphasize that these conclusions are based on the detailed analysis of a single simulation and that further work is needed to put these results on a firmer basis as well as to establish their true applicability to the origin of the dynamical components of the Milky Way. We plan to explore these issues in detail in future papers of this series.

This work has been supported by grants from the U.S. National Aeronautics and Space Administration, the Natural Sciences and Engineering Research Council of Canada, and Fundación Antorchas from Argentina.


format: text

text:
# Optimized Color Gamuts for Tiled Displays

## 1 Introduction

Demand for increased image resolution has grown explosively in recent years, both in imaging devices such as cameras and in displays for those images. In particular, collaborative workspaces require large display walls that can accomodate both presentation-sized graphics and higher resolution images for up-close viewing. In order to satisfy this demand for display resolution, researchers and corporations have resorted to tiled displays, in which multiple display units (typically rear projection systems) each cover a small portion of the overall display area . However, achieving a seamless appearance with such a system poses formidable technical challenges, involving careful geometric alignment of the individual display units, feathering or masking seams between units, correction for distortion, vignetting, and other defects of individual units’ projection systems, and matching of the color characteristics among the different units in the system.

In this paper we address one of these challenges: finding an appropriate color gamut, achievable by all projectors in the system, so that scenes can be rendered on the system without distracting color shifts on the boundaries between different projectors’ output. This problem was previously considered by Stone , who described a non-automated procedure for finding a large common color gamut for tiled displays. Majumder et al.  also studied color matching problems in tiled displays, however their paper addresses only the related problem of mapping a suitable device-independent gamut into individual projectors’ color spaces, and does not discuss how to find such a gamut. We show how to formalize this challenge as a computational geometry optimization problem, and describe efficient algorithms for solving various formulations of this problem.

## 2 Color Gamuts

We review briefly the measurement and algebra of additive color devices such as LCD projectors; for more details see Stone  or her references.

In these devices, pixels are combinations of separate red, green, and blue light signals, each of which has its color determined by filtering light from a white light source, and its intensity determined by a single liquid crystal element. When measured in a suitable device-independent color space, color values can be viewed as elements of a three-dimensional linear algebra: the superposition of two color signals has a value measured by the three-dimensional vector sum of the corresponding two color values. Thus, the gamut of a three-signal additive color device (that is, the set of colors that it can display) naturally forms a parallelepiped in color space, in which the eight parallelepiped corners correspond to eight color values: black (which we denote $`K`$ for short), white ($`W`$), red ($`R`$), green ($`G`$), blue ($`B`$), cyan ($`C`$), magenta ($`M`$), and yellow($`Y`$). However, it is not usually the case that the black point $`K`$ lies at the zero point of the vector space of colors, because a projector displaying black still projects some light onto its screen. Given a suitable set of four of these vectors (say, $`K`$, $`R`$, $`G`$, and $`B`$), we can determine the other four extremal colors as linear combinations of these four: $`W=K+(RK)+(GK)+(BK)=R+G+B\mathrm{2}K`$, $`Y=K+(RK)+(GK)=R+GK`$, etc.

To provide a simplified view of this three-dimensional color space, our figures illustrate two-dimensional analogues. For instance, a person with red-green color blindness can only sense two of the three color bands seen with normal vision; color for such a person would naturally be measured in a two-dimensional device-independent color space, formed by projection from the usual three-dimensional color spaces, and indicated by the coordinate axes in the figures. Similarly, a projector with a short circuit in the green color channel’s signal would generate only a two-dimensional color gamut, in the form of a parallelogram with four corners $`K`$, $`R`$, $`B`$, and $`M=R+BK`$: see Figure 1, left.

We note that some devices, including projector screens based on color wheels and mirror arrays instead of liquid crystals, may form their pixels from the superposition of four or more color signals; for instance certain mirror array displays have a fourth white position on their color wheels . If these signals are independently controllable, the resulting gamut has the form of a zonotope  such as a rhombic dodecahedron; to return to the two-dimensional simplified color model of our illustrations, a viewer with red-green color blindness viewing a three-color projector would see a gamut in the form of a hexagon in his two-dimensional color space (Figure 1, right). In practice, however, the white color is added in progressively as a function of the intensity of the input signal, resulting in gamuts in the form of elongated parallelepipeds. Our algorithms can easily be adapted to handle such shapes without significantly increasing their complexity.

Due to variations in the manufacturing process for their filters, differences in the batch and age of the light bulbs used in their white light sources, and other factors, different projectors, even of the same model, will have gamuts that form different parallelepipeds in the device-independent color space. If one uses these projectors in a tiled display without correcting for their different gamuts, the resulting difference in visual appearance will mar the appearance of seamlessness that most tiled display systems aim for. Thus, we must find a standard gamut that contains only colors that can be displayed by all projectors in our system.

The most obvious approach to this standard gamut problem, from a computational geometry point of view, is to find the three-dimensional shape formed by intersecting the gamuts of all the individual projectors; if we have $`n`$ projectors, this shape can be represented as the intersection of $`\mathrm{6}n`$ halfspaces (the faces of the parallelepipeds from each projector), and constructed in time $`O(nlogn)`$. However, the resulting shape will likely be irregular and not itself be a parallelepiped (Figure 2, left). Non-parallelepiped gamuts commonly arise in subtractive color devices such as inkjet printers, and there has been much work on mapping device-independent colors onto such gamuts , but there are significant technical difficulties in this approach: Complicated algorithms are required to map colors from a device independent space (such as the sRGB  typically used for computer graphics) into such gamuts, and the transformation from the device independent space to the gamut may introduce clipping or color shifts for saturated colors near the boundaries of the gamut. There are difficulties even in using a relatively simple non-parallelepiped gamut such as a warped parallelepiped with curved, bilinear sides; for example, computer graphics algorithms such as anti-aliasing rely on the additivity of colors.

Instead, following Stone , we wish to treat the tiled display system as if it were a single additive color device, by finding a large parallelepiped contained in the intersection of the projector gamuts (Figure 2, right), and using this parallelepiped as our standard gamut. Each projector can represent all the colors in such a standard gamut, and the transformation from the standard gamut in the device independent color space to the gamuts of the individual projectors can be performed in a particularly simple way, by a linear transformation represented by a $`\mathrm{4}\times \mathrm{4}`$ matrix .

Thus, we can state formally the problem of finding a standard gamut as a computational geometry optimization problem: given a suitable objective function $`q`$ (measuring the quality of a given gamut), and a set of $`n`$ parallelepipeds (projector gamuts), find a parallelepiped $`P`$ contained in the intersection of the projector gamuts and maximizing $`q(P)`$. A natural choice for $`q`$ would be the volume of $`P`$, but this seems to lead to inefficient algorithms due to the high number of degrees of freedom, 12, in the specification of $`P`$. Instead it seems appropriate to apply colorimetric principles to simplify the problem as well as to achieve a better match between the objective function and viewer-perceived quality. In the following sections, we explore the choice of objective function and its algorithmic implications.

## 3 Stone’s Procedure

Stone  described a hand-computation procedure for finding large standard gamuts. Her procedure is based on factoring the three-dimensional device-independent color space into a two-dimensional chromaticity and one-dimensional luminosity. The luminosity of a color value is measured by a positive linear combination of the values of its channels; sets of colors with given fixed luminosity values form parallel planes in the device independent color space (Figure 3, left). From the parallelepiped representing any additive color device’s gamut, we can form a chromaticity diagram by translating the black point to the origin of the color space and then projecting the translated gamut onto a plane of constant luminosity; the projection of the parallelepiped forms a triangle in chromaticity space. The inverse image of any point in the projection plane is a line through the origin with constant chromaticity (Figure 3, right).

Stone’s procedure consists of the following steps:

1. Intersect the chromaticity diagrams of the projectors, and find a large triangle within this intersection. Stone does not specify how this triangle is selected, but efficient algorithms are known for finding large triangles within convex polygons .
2. Select the chromaticity of the black point $`K`$ of the standard gamut by averaging the chromaticities of the black points of the individual projectors. Select $`K`$ itself as the darkest value of that chromaticity within the gamuts of all projectors.
3. Select the chromaticity of the white point $`W`$, to be the same as that of the black point and select $`W`$ as the brightest value of that chromaticity within the gamuts of all projectors.
4. Form a standard gamut by using the chromaticity diagram selected in step 1, the black point selected in step 2, and the white point selected in step 3.
5. The resulting gamut may have bright colors that fall outside the gamuts of the individual projectors. If so, scale the gamut by reducing the luminosities of the primary colors until it is contained in the projector gamuts.

Although all of these choices are well motivated, it is difficult to find a quality measure that they optimize, and some of the steps in this procedure would need to be more precisely specified in order to transform it into a computer algorithm.

## 4 Constrained Volume Maximization

Following Stone, we feel that it is appropriate to first select the black and white points of the standard gamut, and then optimize the color. Not only does this help reduce the number of degrees of freedom in the optimization procedure from twelve to six, but it makes sense from a human-computer interaction perspective as well: the prevalence of black and white text in computer displays makes it especially important to achieve a sharp visual distinction between the two colors. To avoid visible color shifts, it is also important that the black and white points have the same chromaticity, as in Stone’s procedure. Thus, we define the following objective function for optimizing a choice of standard gamut:

* Select $`K`$ and $`W`$ to be two points in the intersection of the projector gamuts, both having the same chromaticities, and maximizing the ratio of the luminosities of $`K`$ and $`W`$.
* Choose a standard gamut with maximum volume among all gamuts that have the given values of $`K`$ and $`W`$ and that are contained within all the projector gamuts.

Both of these choices are somewhat arbitrary, and as we will see could be replaced by other similar objectives without significantly changing our optimization algorithms. For instance, other objectives than the ratio of luminosities could be used for selecting $`K`$ and $`W`$, or they could be set independently of each other to minimize and maximize their luminosity values. Compared with Stone’s approach, though, we have the following important differences. First, the initial optimization of $`K`$ and $`W`$ is not altered by a final scaling step; instead, we set $`K`$ and $`W`$ once and for all and then optimize the remaining parameters subject to those constrained values. Secondly, we have an explicit objective function that we seek to optimize (the volume of the gamut), stated independently of any algorithm for performing that optimization.

### 4.1 Setting Black and White

We now describe how to find values of $`K`$ and $`W`$ implementing the first part of our objective function for the optimal standard gamut.

Let $`P`$ denote the polyhedron, in three-dimensional device-independent color space, formed by intersecting the parallelepipeds formed by the color gamuts of the projectors in our tiled display. If we intersect $`P`$ with a line $`\mathrm{}_c`$ through the origin, representing all points with a fixed chromaticity $`c`$, then this line will intersect $`P`$ at two points: $`\lambda ^{}(c)`$ is the closest point to the origin in $`P\mathrm{}_c`$, and has the minimum luminosity of any point in $`P`$ with color $`c`$; $`\lambda ^+(c)`$ is the farthest point from the origin in $`P\mathrm{}_c`$, and has the maximum luminosity of any point in $`P`$ with color $`c`$. Our task is to choose $`c`$ in such a way as to maximize the ratio of the luminosities of $`\lambda ^+(c)`$ and $`\lambda ^{}(c)`$.

To see how to do this, consider projecting the gamut intersection $`P`$ onto a plane of constant luminosity, forming a convex polygon $`\stackrel{\textcolor[rgb]{1,0,0}{bar}}{P}`$ within the triangular chromaticity diagram of our color space. Partition the facets of $`P`$ into two subsets $`P^{}`$ (the facets visible from the origin, containing points $`\lambda ^{}(c)`$) and $`P^+`$ (the facets on the far side of the origin from $`P`$, containing points $`\lambda ^+(c)`$). Then projecting $`P^{}`$ onto the chromaticity diagram forms a subdivision of $`\stackrel{\textcolor[rgb]{1,0,0}{bar}}{P}`$ into smaller convex polygons (Figure 4, left), and projecting $`P^+`$ onto the chromaticity diagram forms a different subdivision of $`\stackrel{\textcolor[rgb]{1,0,0}{bar}}{P}`$ into convex polygons (Figure 4, center). Overlaying the two projections produces a third, more complicated diagram $`\stackrel{\textcolor[rgb]{1,0,0}{hat}}{P}`$ (Figure 4, right).

###### Lemma 1

The optimal value of $`c`$ is one of the vertices of diagram $`\stackrel{\textcolor[rgb]{1,0,0}{ℎ𝑎𝑡}}{P}`$.

The proof idea is that within any of the higher dimensional faces of the diagram, the value we are trying to optimize is a fractional linear function that cannot have any local optima; we omit the details. Thus, we are led to the following algorithm for finding $`K`$ and $`W`$, which implicitly enumerates all the vertices of this diagram.

1. Construct a representation of the vertices, edges, and facets of polyhedron $`P`$ by intersecting the gamuts of the individual projectors. Partition the boundary of $`P`$ into $`P^{}`$ and $`P^+`$.
2. For each vertex $`v`$ of $`P`$ in $`P^{}`$, let $`c_v`$ be the chromaticity of the vertex; the vertex itself is $`v=\lambda ^{}(c_v)`$. Construct the line through the origin with fixed chromaticity $`c_v`$, and intersect this line with the planes through each facet of $`P^+`$. Let $`\lambda ^+(c_v)`$ be the nearest of these intersection points to the origin.
3. Similarly, for each vertex $`w=\lambda ^+(c_w)`$ in $`P^+`$, construct the line with fixed chromaticity $`c_w`$, intersect this line with the planes through each facet of $`P^{}`$, and let $`\lambda ^{}(c_w)`$ be the farthest of the resulting intersection points from the origin.
4. For each two line segments $`sP^{}`$ and $`tP^+`$, project these two segments into the chromaticity diagram. If they intersect in a single point $`c_{st}`$, let $`\lambda ^{}(c_{st})`$ denote the point of chromaticity $`c`$ on $`s`$, and let $`\lambda ^+(c_{st})`$ denote the point of chromaticity $`c`$ on $`t`$.
5. Compare the luminosity ratios of all the pairs $`\{\lambda ^{}(c),\lambda ^+(c)\}`$ constructed in steps 2–4, let $`c^{}`$ denote the value of $`c`$ maximizing this ratio among these pairs, let $`K=\lambda ^{}(c^{})`$, and let $`W=\lambda ^+(c^{})`$.

###### Lemma 2

The algorithm described above finds the optimal values of $`K`$ and $`W`$ in $`O(n^\mathrm{2})`$ time.

We note that a similar diagram overlay technique can be used to find the optimal choice of $`K`$ and $`W`$ (constrained to have a common chromaticity) for many other objective functions than our choice of the luminosity ratio. It seems likely that the time for these steps can be reduced using techniques for optimization in overlays of two planar diagrams , however this would not reduce the complexity of our overall algorithm.

### 4.2 Filling out the gamut

We now discuss the more difficult algorithmic problem of computing the maximum-volume gamut for a fixed choice of $`K`$ and $`W`$. As discussed in Section 2, this gamut can be determined by a suitable set of four of its eight extremal colors. Once $`K`$ and $`W`$ are chosen, we need only fix two more noncomplementary colors, say $`R`$ and $`B`$; all other extremal colors can be computed by formulas that are linear in the coordinates of $`R`$, $`K`$, $`B`$, and $`W`$. However, $`R`$ and $`B`$ can be chosen independently among the points of the device-independent color space. Thus, the set of possible gamuts that we must search forms a six-dimensional space in which $`R`$ and $`B`$ each contribute three coordinates. We denote points in this space by pairs $`(R,B)`$ where each of $`R`$ and $`B`$ are themselves three-dimensional colors.

###### Lemma 3

In the six-dimensional space $`(R,B)`$, the set of gamuts that are contained in the intersection of all projector gamuts corresponds to a polytope $`\mathrm{\Gamma }`$ formed by the intersection of $`\mathrm{3}\mathrm{6}n`$ halfspaces, where $`n`$ denotes the number of projectors in the system.

###### Proof

If $`P`$ and $`Q`$ are any two convex polytopes, then $`PQ`$ if and only if the vertices of $`P`$ all belong to $`Q`$. Thus, a standard gamut determined by the four extremal colors $`K`$, $`W`$, $`R`$, and $`B`$ is contained within the intersection of the printer gamuts if and only if all eight extremal colors are contained within this intersection. The two colors $`K`$ and $`W`$ are guaranteed by our contruction to already be within this intersection. For any other color $`c\{R,G,B,C,M,Y\}`$, any projector gamut $`\mathrm{\Gamma }_i`$, and any facet $`f`$ among the six facets of $`\mathrm{\Gamma }_i`$, the plane through $`f`$ must not separate $`c`$ from $`\mathrm{\Gamma }_i`$. Since each color $`c`$ can be expressed as a linear combination of the coordinates of $`R`$ and $`B`$ (and of the constant coordinates of $`K`$ and $`W`$), the constraint that $`c`$ and $`\mathrm{\Gamma }_i`$ must be on the same side of $`f`$ can be expressed as a linear inequality in the six-dimensional space $`(R,B)`$. There are six choices of $`c`$, $`n`$ choices of $`\mathrm{\Gamma }_i`$, and six choices of $`f`$ for each $`\mathrm{\Gamma }_i`$, so the total number of these linear inequalities is $`\mathrm{3}\mathrm{6}n`$.

We next consider the nature of the objective function we are trying to optimize. On the face of it, this is a cubic function: the determinant of the $`\mathrm{3}\times \mathrm{3}`$ matrix the rows of which are the three-dimensional vectors $`RK`$, $`BK`$, and $`GK`$. However, we can simplify this considerably: by additivity of determinants, and the fact that $`GK=(WK)(RK)(BK)`$, we can separate the volume determinant into a sum of three determinants of simpler matrices $`vol(R,B)=det(RK,BK,WK)det(RK,BK,RK)det(RK,BK,BK)=det(RK,BK,W+K)`$. Since $`WK`$ is already fixed, the function to maximize over $`\mathrm{\Gamma }`$ is therefore seen to be a quadratic over the coordinates $`(R,B)`$. Unfortunately, this quadratic is neither convex nor concave. If it were convex, we could find the optimal standard gamut as a convex program, using generalized linear programming techniques, in time $`O(n)`$ , while if it were concave, the optimal gamut would be a vertex of $`\mathrm{\Gamma }`$, and could be found by enumerating all vertices of $`\mathrm{\Gamma }`$ and computing the objective function for each.

However, the quadratic nature of $`vol(R,B)`$ does help us somewhat: suppose we want to determine the maximum of $`vol(R,B)`$ within a compact convex set $`S`$. Find the affine hull of $`S`$ (that is, the smallest linear subspace of the six-dimensional space $`(R,B)`$ that contains all of $`S`$), and restrict the function $`vol`$ to this subspace. If the restricted quadratic is positive definite or indefinite, it is unbounded in the linear subspace, and the maximum must occur on the boundary of $`S`$. If the restricted quadratic is negative definite, it has a unique local maximum on the linear subspace. If this maximum occurs within the interior of $`S`$, it is the maximum on $`S`$ as well as on the larger linear subspace. Otherwise, the maximum again must occur on the boundary of $`S`$. We can use the facts in the following algorithm for finding $`(R,B)`$ maximizing $`vol(R,B)`$:

1. Find the set of $`\mathrm{3}\mathrm{6}n`$ 6-dimensional halfspaces, described in Lemma 3, such that any pair $`(R,B)`$ contained within the intersection of these halfspaces determines a feasible standard gamut.
2. Construct the polytope $`\mathrm{\Gamma }`$ formed by intersecting the halfspaces. The polytope’s worst case complexity is $`O(n^\mathrm{3})`$ by McMullen’s upper bound theorem , and worst-case optimal halfspace intersection algorithms (by projective duality, equivalent to convex hull construction) with time matching this complexity are known .
3. Subdivide the boundary of $`\mathrm{\Gamma }`$ into simplices (triangles, tetrahedra, and their higher dimensional analogues). This can be done by recursively subdividing all lower dimensional boundary faces, and then within each facet connecting the bottommost vertex to each simplex in the triangulation of the facet’s boundary. This pulling triangulation again has $`O(n^\mathrm{3})`$ complexity by the upper bound theorem.
4. Within each $`k`$-dimensional simplex $`\mathrm{\Delta }`$ of the triangulation, restrict the function $`vol(R,B)`$ to the affine hull of $`\mathrm{\Delta }`$ by substituting linear equations for $`\mathrm{6}k`$ of the coordinates of $`(R,B)`$. Test whether the resulting restricted quadratic function is negative definite. If so, and the maximum of the function is within $`\mathrm{\Delta }`$, then let $`(R_\mathrm{\Delta },B_\mathrm{\Delta })`$ denote the point achieving this maximum.
5. Compute $`vol(R_\mathrm{\Delta },B_\mathrm{\Delta })`$ for all points $`(R_\mathrm{\Delta },B_\mathrm{\Delta })`$ constructed in the previous step, and let $`(R,B)`$ denote the pair achieving the maximum volume.

Other objective functions could be substituted, as long as we can find optimal values of the function within six-dimensional simplices, without increasing the overall complexity of the algorithm.

### 4.3 Summary of Algorithm

To summarize, our volume maximization algorithm consists of finding $`K`$ and $`W`$ by overlaying two chromaticity diagrams as described in Section 4.1, finding $`R`$ and $`B`$ by searching for maxima of $`vol(R,B)`$ within each simplex of a triangulation of the six-dimensional polytope $`\mathrm{\Gamma }`$ as described in Section 4.2, and using $`K`$, $`W`$, $`R`$, and $`B`$ to determine the standard gamut.

###### Theorem 4.1

The algorithm described above finds a standard gamut optimizing our objective function in time $`O(n^\mathrm{3})`$, where $`n`$ denotes the number of projector gamuts in the input.

We note that further speedups are possible if we allow for approximate results: a $`(\mathrm{1}ϵ)`$-approximation to the volume-maximizing gamut for given corners $`K`$ and $`W`$ can be found in time $`O(nlogn+ϵ^\mathrm{3})`$ by replacing the intersection of the projector gamuts with an approximate intersection having $`O(\mathrm{1}/ϵ)`$ vertices. This approximate intersection can be found by using an affine transformation to ensure that the true intersection is sufficiently fat, then applying an algorithm of Dudley  to replace the intersection by a polytope with fewer vertices; we omit the details.

## 5 Quasiconvex Programming

Our discussion above centers on the optimization of objective functions involving volume of the color gamut. However, our algorithms could easily be adapted to other objective functions, and it is not clear that volume is the correct choice. We see two potential problems with our volume-based approach. First, the running time, while cubic, has a high constant factor due to the factor of $`\mathrm{3}\mathrm{6}`$ in Lemma 3; this may make our running times impractically slow. This speed issue could be addressed from the point of view of approximation algorithms, as in the previous section, but it might be more appropriate to find a faster exact algorithm for an easier-to-optimize quality measure. Second, volume maximization neglects another issue of color gamut optimization: how do we match the corners of the chosen standard gamut to the usual additive primary colors $`R`$, $`G`$, and $`B`$? Any permutation of these corners would produce a geometrically equivalent parallelepiped but only one of these six permutations is likely to be appropriate for use as our standard gamut. In practice, we expect the choice of corner labeling to be obvious, but this ambiguity indicates a possible weakness of our volume-based approach. In this section, we present a more generic approach based on quasiconvex programming that would have the advantages of faster running time and automatic gamut corner labeling; however more colorimetric expertise would be needed to make specific design choices among those offered by this approach.

Quasiconvex programming is a general framework for LP-type computational geometry optimization problems, formulated by Amenta et al.  in the context of mesh smoothing procedures, and later used by the authors for several applications in information visualization and mesh generation . Define a nested convex family to be a function $`\kappa (t)`$ mapping real values $`t`$ to convex sets in some Euclidean space $`E^d`$, such that for any $`t<t^{}`$, $`\kappa (t)\kappa (t^{})`$, and such that for all $`t`$, $`\kappa (t)=_{t^{}>t}\kappa (t^{})`$. If $`\kappa `$ is a nested convex family, we can define a quasiconvex function $`f_\kappa (x)=inf\{t|x\kappa (t)\}`$ mapping $`E^d`$ back to the real numbers; conversely every quasiconvex function corresponds to a nested convex family. An instance of a quasiconvex program consists of a finite set $`K=\{\kappa _\mathrm{0},\kappa _\mathrm{1},\mathrm{}\}`$ of nested convex families (or equivalently quasiconvex functions). The desired output is the pair

$$inf\left\{(t,x)|x\underset{\kappa _iK}{}\kappa _i(t)\right\}$$

where the infimum is taken in the lexicographic ordering, first by $`t`$ and then by the coordinates of $`x`$. Less formally, we seek the point $`x`$ minimizing the maximum value of $`f_{\kappa _i}(x)`$.

As Amenta et al.  showed, when $`d`$ is a fixed constant, any quasiconvex program consisting of $`n`$ nested convex families can be solved by generalized linear programming techniques, assuming the existence of primitives for evaluating the input quasiconvex functions as well as for solving subproblems with a constant number of inputs. The time for these techniques is $`O(En+Slogn)`$ where $`E`$ is the time for an evaluation and $`S`$ is the time for a subproblem solution. Alternatively, Amenta et al. showed that simple hill-climbing techniques for quasiconvex programs are guaranteed to converge to the globally optimal value.

To formulate the color gamut optimization problem as a quasiconvex program, recall that additive gamuts can be parametrized by a four-tuple of colors, say $`(K,R,G,B)`$, and that the other four colors $`C,M,Y,W`$ can be expressed as linear functions of these four colors. We need functions $`q_c(x)`$ for each gamut corner $`c\{K,R,G,B,C,M,Y,W\}`$, measuring how well a point $`x`$ in device independent color space matches some desired ideal location for $`c`$. These functions $`q_c`$ could be Euclidean distances from ideal locations, or linear objective functions measuring distance along ideal directions (”blackest black”, ”reddest red”, and so forth). It may be appropriate to assign greater weight to the locations of $`K`$ and $`W`$ than to the other corners. Lower values of $`q_c`$ will indicate a better fit to the ideal location, and higher values will indicate a worse fit. We will use these functions to measure the quality of a gamut $`\mathrm{\Gamma }`$, in the 12-dimensional space of color gamuts parametrized by color 4-tuples $`(K,R,G,B)`$, where the quality $`Q_c(\mathrm{\Gamma })`$ is measured by finding the point labeled $`c`$ in gamut $`\mathrm{\Gamma }`$ and applying $`q_c`$ to that point.

###### Lemma 4

$`Q_c`$ is quasiconvex if and only if $`q_c`$ is quasiconvex.

Thus, in order to apply quasiconvex programming, all we need is for each of the functions $`q_c`$ to be quasiconvex.

###### Theorem 5.1

Suppose that for each $`c\{K,R,G,B,C,M,Y,W\}`$ we are given a quasiconvex function $`q_c`$ on device-independent color space, and define $`Q_c`$ as above. Then the problem of finding the gamut $`\mathrm{\Gamma }`$ minimizing $`max_cQ_c(\mathrm{\Gamma })`$ can be solved in linear time by quasiconvex programming.

It is natural to hope that colorimetric principles will lead to a good choice of the exact form of the functions $`q_c`$; we leave the details of this choice for further work. We note that the formulation above does not include any constraint that $`W`$ and $`K`$ have the same chromaticities; it would be of interest to incorporate such a constraint into our quasiconvex programming framework, or to determine whether the color shifting it prevents would actually be a problem in the unconstrained quasiconvex program.

## 6 Conclusions and Further Work

We have formulated problems of color gamut intersection as that of optimizing an objective function on parallelepipeds within a convex polyhedron, and found efficient algorithms for this optimization task for several colorimetrically motivated definitions of the objective function. The algorithms described in this paper have not yet been implemented, but we believe they should be simple and efficient enough to be suitable for typical problem sizes; we hope to continue this work by implementing some of our algorithms and comparing their results with those of Stone’s procedure.

More generally, color spaces seem likely to provide a fruitful source of problems in three-dimensional computational geometry. Other problems that have been studied in this area include color gamut mapping , color quantization , and optimized printing using restricted numbers of ink colors  but it seems likely that much more work remains to be done in the computational geometry of color.

## Acknowledgements

Work of Eppstein was done in part while visiting Xerox PARC and supported in part by NSF grant CCR-9912338. We thank Maureen Stone for helping us understand the mysterious world of color.


text:
# 1 Current for (a) a quantum dot and (b) a junction with 𝑉=0, defined by Eqs. (2) and (5) in Ref. , , respectively. The parameters are the same as in Ref. , .

Comment on “Time-Dependent Density-Matrix Renormalization Group: A Systematic Method for the Study of Quantum Many-Body Out-of-Equilibrium Systems”

In a recent Lettercm , Cazalilla and Marston (CM) proposed a time-dependent density-matrix renormalization group (TdDMRG) algorithm for accurate evaluation of out-of-equilibrium properties of quantum many-body systems. An interesting result they obtained is that for a point contact junction between two Luttinger liquids, a current oscillation develops after initial transient in the insulating regime (Fig. 2 in Ref. 1). They attributed this oscillation to a non-linear response of the Luttinger liquid. Here we would like to point out that (a) the oscillation they observed is an artifact of the method; (b) the TdDMRG can be significantly improved by extending the definition of the density matrix to adapt the non-equilibrium evolution of the ground state.

In the TdDMRG scheme of CM, it is the equilibrium ground state that is targeted and used in the construction of the reduced density matrix. At the beginning of the evolution, the state does not deviate much from the ground state $`|\mathrm{\Psi }_0`$, and results obtained with that approach are accurate. However, the time evolution of the ground state evokes excited states in a non-equilibrium system, so the long time behavior of the state becomes very poor and could not be substantially improved by keeping more states. This can be clearly seen from Fig. 1, where the numerical results obtained with that approach (i.e. the $`N_t=0`$ curves) for the two models in Ref. cm,  begin to deviate from the exact ones at $`t25`$.

To retain the information on the relevant excitation states, we propose to define the density matrix from the time dependent wavefunction $`|\mathrm{\Psi }(t)`$ as

$$\rho =\underset{i=0}{\overset{N_t}{}}\alpha _i|\mathrm{\Psi }(t_i)\mathrm{\Psi }(t_i)|,\underset{i}{}\alpha _i=1,$$
(1)
where$`|\mathrm{\Psi }(t_0)=|\mathrm{\Psi }_0`$, $`t_0`$ is the starting time at which a bias is imposed and $`N_t`$ is the number of equal time intervals within the whole evolution time to be evaluated. The reduced density matrix for the system is defined by tracing over the degrees of freedom of the environment. Both the time evolution of the wavefunction and the newly defined reduced density matrix should be evaluated from the beginning of the DMRG truncation, not just for the final lattice system. The density matrix used by CM corresponds to the $`N_t=0`$ case. In our calculation, we took $`\alpha _i=1/2`$ for $`i=0`$ and $`1/2N_t`$ otherwise. Empirically, we find that $`N_t6`$ is generally needed. Fig. 1 shows our results with $`N_t=24`$. By comparison with the exact results, we find that this new scheme significantly improves the large time scale behavior of the current.

Now let us turn to the tunnelling junction between two Luttinger liquids in the insulating regime. For $`N_t=0`$, our calculation confirms the existence of a current oscillation around $`t20`$. However, we find that the oscillation depends on the number of states retained. By further applying the new TdDMRG scheme to this system, we find that the oscillation does not exist at all. In Fig. 2, the numerical results with $`N_t=0`$ and $`N_t=24`$ are shown for detailed comparisons. Therefore, the current oscillation observed by CM is in fact an artifact of the TdDMRGcm , rather than a nonlinear response of the system to the applied bias. This also shows that in the study of physical properties of a quantum state out of equilibrium, it is important to include the relevant excitation states in the definition of the reduced density matrix.

PACS number: 71.27.+a, 71.10.Pm, 72.15.Qm, 73.63.Kv

H.G. Luo1,2,T. Xiang1, and X.Q. Wang1

1Institute of Theoretical Physics
Chinese Academy of Sciences
P.O.Box 2735, Beijing 100080, China

2Department of Modern Physics
Lanzhou University
Lanzhou 730000, China