Proceedings ofthe 29th Grain Formation Workshop /Dust in Galaxies 2011

held at Center of Planetary Science (CPS), Kobe, Japan
9—11 November 2011

Eds.: S. Okuzumi, I. Sakon, Y. Kimura, A. Inoue, H. Kimura, H. Miura, T. Muto, T. Yamamoto
Online materials; https://www.cps-jp.org/~gfw/pub/2011-11-09/

1Proceedings ofthe 29th Grain Formation Workshop /Dust in Galaxies 2011

held at Center of Planetary Science (CPS), Kobe, Japan
9—11 November 2011

Eds.: S. Okuzumi, I. Sakon, Y. Kimura, A. Inoue, H. Kimura, H. Miura, T. Muto, T. Yamamoto
Online materials; https://www.cps-jp.org/~gfw/pub/2011-11-09/

2

Program
9th November, 2011
13:30-13:40 Tetsuo Yamamoto (Hokkaido University) / Akio Inoue (Osaka Sangyo University)“Opening Remarks”
13:40-14:00 Akemi Tamanai (Heidelberg University)“A study of dust grain formation: A role of minor elements”
14:00-14:20 Chiyoe Koike (Osaka University)“IR spectra of silica (SiO2) polymorphs”
14:20-14:40 Steffen Wetzel (Heidelberg University)“Condensation and solid phase reactions of Mg silicate system”(coffee break)
15:00-15:40 Aki Takigawa (University of Tokyo)“Morphological analyses of presolar alumina grains towards understanding the circumstellar
dust formation and evolution processes”
15:40-16:00 Yuki Kimura (Tohoku University)“Homogeneous nucleation experiment of cosmic dust”(coffee break)
16:20-16:40 Akane Sakurai (Nagoya university)“Star formation and dust extinction properties of local galaxies seen from AKARI”
16:40-17:00 Tsutomu Takeuchi (Nagoya University)“A Study on the Star Formation History of Galaxies Revealed by the UV-IR Bivariate Luminosity
Function”(coffee break)
17:20-17:40 Shogo Masaki (Nagoya University)“An analytic model for the distribution of matter and dust around galaxies”
17:40-18:00 Daisuke Yamasawa (Hokkaido University)“From dust to galaxy in the early universe”
18:00-18:20 Ryosuke Asano (Nagoya University)“Evolution of grain size distribution of galaxies”

3

10th November, 2011
9:20-9:40 Hitoshi Miura (Tohoku University)“A new diagnostic method for growth condition of rock-forming minerals”
9:40-10:00 Nagisa Machii (kobe University Graduate school of science)“Conditions for chondrules to be intruded and embedded into matrix”(coffee break)
10:20-10:40 Masafumi Matsumura (Kagawa University)“Verification on the Alignment Mechanism of Dust Grains by Radiative Torques:Correlation between the Interstellar Polarization and Dust Temperature”
10:40-11:00 Mitsuyoshi Yamagishi (Nagoya University)“An analytic model for the distribution of matter and dust around galaxies”
11:00-11:20 Mina Ishikura (Hokkaido University)“Lyman alpha transfer in multiphase system - effect of internal motion of cloud”-- poster presentation (5min talk) -11:20-11:25 Chihiro Kaito (Ritsumeikan University)“Laboratory analogy of amorphous enstatite fine grain formation and crystallization”
11:25-11:30 Takuya Ikeda (Kyoto University)“Development of Mercury Dust Monitor (MDM) for BepiColombo MMO”
11:30-11:35 Itsuki Sakon (University of Tokyo)“Mid-Infrared Multi-Epoch Observations of Dust Forming Nova V1280Sco”
11:35-11:40 Okuzumi Satoshi (Nagoya University)“Coagulation and Radial Drift of Porous Dust Aggregates in Protoplanetary Disks”(lunch break)
13:20-14:00 Akira Tsuchiyama (Osaka University)“Three-dimensional structure of Itokawa regolith particles.”
14:00-14:40 Eizo Nakamura (Okayama University)“Comprehensive analyses of micro-grains from the surface of asteroid returned by Hayabusa:
Implications for collisional space environments”(coffee break)
15:10-15:30 Hiroki Senshu (Chiba Institute of Technology Planetary Exploration Research Center)“Levitation dust: a source of IDPs”
15:30-15:50 Junya Matsuno (Osaka University)“An synthesis experiment of GEMS analogue material produced by thermal plasma”

4(coffee break)
16:10-16:50 Shota Nunomura (National Institute of Advanced Industrial Science and Technology)“Coagulation and growth of Dust in Plasma by Charge Fluctuations ~from Molecules to clusters
and dust”
16:50-17:10 Yasuaki Hayashi (Kyoto Institute of Technology Faculty of Engineering and Design)“Coulomb-Crystallization and Behavior of Fine Particles in Experimental Plasmas”
11th November, 2011
9:20-9:40 Koji Wada (Chiba Institute of Technology Planetary Exploration Research Center)“Numerical simulation of collisions between dust aggregates consisting of particles with a size
distribution”
9:40-10:00 Hiroshi Kobayashi (Friedrich Schiller University Jena)“Necessary condition for gas giant formation”
10:00-10:20 Yuri Fujii (Nagoya University)“Ionization Degrees with Dust Grains in Circumplanetary Disks”(coffee break)
10:40-11:20 Masashi Arakawa (Kyushu University)“Hydrogen ordering in ice observed from neutron diffraction: Application to planetary science”
11:20-11:40 Kyoko Tanaka (Hokkaido University, Institute of Low Temperature Science)“Evaporation of icy planetesimals due to planetesimal bow shocks”(lunch break)
13:40-14:20 Atsuo Okazaki (Hokkai-Gakuen University)“Structure of the Wind-Wind Collision Shocks in Massive Binaries”
14:20-14:40 Toshiya Ueta (University of Denver Division of Natural Science and Mathematics)“Interactions between stellar winds and the interstellar medium and their effects on dust grains”(coffee break)
15:00-15:20 Toru Kondo (Nagoya University)“Properties of dust and PAHs in clouds with high-energetic phenomena based on the AKARI
all-sky-survey maps”
15:20-15:40 Kohji Tsumura (ISAS/JAXA)“Distribution of the 3.3µm PAH band in our Galaxy detected by the AKARI IRC”

5

Abstract“A study of dust grain formation: A role of minor elements”
TAMANAI, Akemi (Heidelberg University)
Large amount of silicate dust grains has been detected in many astronomical environments.
However, silicate dust grain formation itself is not well understood yet. Although theoreticalsimulations for possible condensates illustrate how abundant the different silicate dust materials arein astronomical surroundings and they give us a key to clarify the formation of silicates, thesesimulations do not cover all dust species, which possibly exist. In order to deepen the understanding
from the first species condensing out of the gas to the entire condensation sequence of dust
materials present in the environments, first we take into account the species, which contain minor
elements (Al, Na, and Li) and perform the mid-IR spectroscopic measurements for obtaining
benchmark spectra.

“IR spectra of silica (SiO2) polymorphs”
KOIKE, Chiyoe (Osaka University)
As for circumstellar dust, Mg-rich silicates such as forsterite and enstatite were mainly detected.
Recently, silica SiO2 is detected in debris disks. This referred silica is various crystal types such as
quartz, fused quartz, annealed silica. We report on the infrared spectra of various types' silica anddiscuss on possibility of silica and observation.

“Condensation and solid phase reactions of Mg silicate system”
WETZEL, Steffen (Heidelberg University)
We have been performing in-situ IR spectroscopy during formation process of layers byco-evaporation of Fe, Mg, MgO, and SiO under well-defined pressure and temperature in ultra highvacuum (UHV) system in order to deepen the understanding of IR spectral changes uponcondensation conditions, especially effects related to elemental Fe particles enclosed in silicates.

6“Morphological analyses of presolar alumina grains towards understanding the circumstellar dust
formation and evolution processes”
TAKIGAWA, Aki (University of Tokyo)
Presolar grains are the condensates around evolved stars, which have survived the isotopic
exchange processes in the early solar system. Alumina (Al2O3), one of the first condensates aroundO-rich evolved stars, is observed by infrared spectroscopy for many AGB stars and present inchondrites as presolar grains. In order to understand the mineralogical nature of presolar alumina
and the formation and evolution processes of circumstellar alumina grains, we performed systematic
analyses of morphology, crystal structure, and oxygen isotopic compositions of alumina grains inunequilibrium ordinary chondrites, and dissolution experiments of alumina polymorphs“Homogeneous nucleation experiment of cosmic dust”
KIMURA, Yuki (Tohoku University)
Homogeneous condensation is a first occurring process to obtain solid materials around evolvedstars. Since there is no or very few heterogeneous nucleation site, solid materials must condense bythemselves homogeneously.

However, there is no data that how large super-saturation is required

for the condensation. To investigate the homogeneous nucleation and growth process of cosmic dust,
interferometric observation was attempted for the fist time to the gas evaporation method inlaboratory.

In preliminary experiment, we found the degree of super-saturation for nucleation was

extremely high, 106-107, which was determined from the interferogram. Surface free energy, sticking
probability and growth velocity can also be deduced from the calculated values by thesemi-phenomenological nucleation theory.

“Star formation and dust extinction properties of local galaxies seen from AKARI”
SAKURAI, Akane (Nagoya university)
Accurate estimation of the star formation -related properties of galaxies is crucial for understanding
the evolution of galaxies. In galaxies, ultraviolet (UV) light emitted by formed massive stars is
attenuated by the dust which is also produced by the SF activity, and is reemitted at mid- and farinfrared wavelengths (IR). In this study, we investigated the star formation rate (SFR) and dust
extinction using data at UV and IR. We selected local galaxies, which are detected at AKARI FIS
90um band. We measured flux densities at FUV and NUV from the GALEX images. We examined
the SF and extinction by using 4 wave bands given by AKARI. Then, we calculated FUV and total IR
luminosities, and obtained the so called SF luminosity (the total luminosity related to star formation
activity) and the SFR. We found that in most of galaxies, SF luminosity is dominated by luminosity of
dust. We also found that galaxies with higher SF activity have a higher fraction of SF hidden by dust.
Especially, SF of galaxies which have SFRs>20Myr is almost completely hidden by dust. Although

7

these results were claimed by previous studies, confirming them precisely using much larger
samples from AKARI and GALEX all sky surveys has a great impact on our understanding of the SF
in Local galaxies.

“A Study on the Star Formation History of Galaxies Revealed by the UV-IR Bivariate Luminosity
Function”
TAKEUCHI, Tsutomu (Nagoya University)
Star formation in galaxies is tightly connected to the dust grain formation through the production ofheavy elements. Observationally, it is the most direct to measure the star formation rate by theamount of ionizing UV photons, but dust formed by the star formation scatter and absorb thesephotons, and finally re-emit the energy at the IR. Hence, it is necessary to observe galaxies at both
UV and IR to explore their star formation in the cosmic history. In this presentation, first we introduce
a statistical method to handle the two datasets, which have been treated independently at eachwavelength in a unified manner, and reveal the star formation in galaxies hidden by dust.

“An analytic model for the distribution of matter and dust around galaxies”
MASAKI, Shogo (Nagoya University)
We develop a simple analytic model to make a quick interpretation of the observed mass and dust
surface density profiles around galaxies reported by Menard et al. (2010). Our model is based on theso-called halo model. We show our model reproduces the observed profiles very well, and discuss
how far dust spreads from galaxies.

“From dust to galaxy in the early universe”
YAMASAWA, Daisuke (Hokkaido University)
We present a semi-analytic model for the formation and evolution of high-redshift galaxies. We
investigate cosmic star formation rates and the reionization history, using our galaxy model in which
we include a novel implementation of dust size evolution and resulting molecular hydrogen (H2)
formation on dust grains in the early universe. We show cosmic star formation efficiency of our model
is low, because the H2 formation rate on dust grains is suppressed by the dust destruction by reverse
shocks in supernova remnants. Our results show need of additional reioniation sources or top-heavy
IMF in z>6 for the cosmological reionization. We conclude that not only the amount but also the size
distribution of dust strongly affects the cosmic star formation efficiency and the IMF transition fromPop III to Pop II.

8“Evolution of grain size distribution of galaxies”
ASANO, Ryosuke (Nagoya University)
To understand the galaxy evolution, it is essential to understand the evolution of dust in galaxies. Up
to now, dust evolution has been studied with various dust evolution models, but most of them arebased on the assumption of a typical size of grains. However, since grain size distribution of galaxies
varies by various processes (dust destruction by SN shocks, metal accretion onto surface of grains
and etc.) [Yamasawa et al. (2011), Hirashita et al. (2011)], grain size distribution is a very important
ingredient to understand the dust evolution. In this talk, we would introduce the dust evolution modelwith grain size distribution and report our first results.

“A new diagnostic method for growth condition of rock-forming minerals”
MIURA, Hitoshi (Tohoku University)
Chondritic meteorites contain rock-forming minerals in which the chemical compositions are notuniform. The in-homogeneity is termed as zoning profile, which is considered to reflect the growth
condition (cooling rate, closed- or open-system, super-cooling, and so forth) when they were formedin early solar system. We investigated the zoning profile, growth rate, and mineral morphology
formed under a constant cooling rate based on a numerical scheme. We derived a new analytic
expression to predict the growth rate and zoning profile with respect to the cooling rate. The analytic
formula can be used as a diagnostic tool to evaluate the growth condition of solid-solution minerals,(e.g., olivine) from its zoning profile.

“Conditions for chondrules to be intruded and embedded into matrix”
MACHII, Nagisa (Kobe University)
Collision experiments were performed using glass bead of 3 mm in diameter as a chondrule analog
and polydisperse silica particles of about 0.8 µm in diameter as a matrix analog. We investigatedimpact velocity and matrix size needed to capture the colliding chondrules. We showed that toembed a chondrule into matrix, impact pressure should be larger than the static strength of matrix.
The minimum size for capturing a chondrule is at least 3 cm in diameter.

“Verification on the Alignment Mechanism of Dust Grains by Radiative Torques: Correlation between
the Interstellar Polarization and Dust Temperature”
MATSUURA, Masafumi (Kagawa University)
N/A9“An analytic model for the distribution of matter and dust around galaxies”
YAMAGISHI, Mitsuyoshi (Nagoya University)
We develop a simple analytic model to make a quick interpretation of the observed mass and dust
surface density profiles around galaxies reported by Menard et al.(2010). Our model is based on theso-called halo model. We show our model reproduces the observed profiles very well, and discuss
how far dust spreads from galaxies.

“Lyman alpha transfer in multiphase system - effect of internal motion of cloud”
ISHIKURA, Mina (Hokkaido University)
For study of high-redshift galaxies, Lyα observation is crucial. Strong Lyα emission observed fromhigh-redshift galaxies is problem, because of large optical depth for Lyα. Neufeld(1991) shows thatassumption of multiphase structure (cloud and inter-cloud medium) permits higher escape fraction.
However he neglects effect of internal motion of cloud. I will talk about Neufeld model, and results ofmy simulation with assumption of internal motion of dusty cloud.

“Laboratory analogy of amorphous enstatite fine grain formation and crystallization”
KAITO, Chihiro (Ritsumeikan University)
Amorphous enstatite (MgSiO3) grains were produced by the simultaneous evaporation of SiO andMg vapor by up and down double heaters method in Ar gas pressure of 10 Torr. Produced particles
were mixture of MgO crystallites and amorphous MgSiO3 structure. High-resolution electron
microscopy showed the crystallites of MgO and MgSiO3 crystallites less than 10nm. Crystallization of
MgSiO3 amorphous structure have been determined by heating in vacuum.

“Development of Mercury Dust Measurement (MDM) for BepiColombo MMO”
IKEDA, Takuya (Kyoto University)
N/A

10

“Mid-Infrared Multi-Epoch Observations of Dust Forming Nova V1280Sco”
SAKON, Itsuki (University of Tokyo)
V1280Sco is a slow dust-forming nova that appeared only at a distance of 1.6kpc. The evidence ofdust formation was reported on 23 days after the discovery. We have carried out mid-infrared
photometric and spectroscopic observations of this nova on day 150 with Subaru/COMICS and ondays 1272 and 1616 with Gemini-S/TReCS. Then the near- to mid-infrared spectral energy
distributions (SEDs) of V1280Sco on days 150, 1272 and 1616 are fitted by multi-components dust
emission model. We have found that the astronomical silicate emission features only appear in the
spectra obtained on 1272 and 1616, suggesting the later formation of silicate grains in the case of
V1280Sco. In this presentation I will discuss the chemical evolution of dust grains formed in theejecta of nova V1280Sco.

“Coagulation and Radial Drift of Porous Dust Aggregates in Protoplanetary Disks”
OKUZUMI, Satoshi (Nagoya University)
Rapid orbital drift of macroscopic dust particles is one of the major obstacles against planetesimalformation in proto-planetary disks. We reexamine this problem by considering porosity evolution ofdust aggregates. We show that dust particles can evolve into highly porous aggregates (with internal
densities of much less than 0.1 g/cc) even if collisional compression is taken into account.
Furthermore, we find that the high porosity triggers significant acceleration in collisional growth.
Thanks to this rapid growth, the highly porous aggregates are found to overcome the radial drift
barrier at distances less than 10 AU from the central star. This suggests that, if collisional
fragmentation is truly insignificant, formation of icy planetesimals is possible via direct collisional
growth of submicron-sized icy particles.

“Three-dimensional structure of Itokawa regolith particles.”
TSUCHIYAMA, Akira (Osaka University)
Particles recovered by the Hayabusa spacecraft are considered to be regolith particles on asteroid
Itokawa. Their 3D structures and external shapes were obtained using x-ray micro-tomography at
SPring-8. The modal abundances of minerals and density of the bulk sample are similar to those of
LL chondrites. Any large-scale melting textures, such as agglutinates in lunar regolith, were notobserved. The size distribution and shape distribution are also different from those of lunar regolith.
Based on these data, formation and evolution of regolith particles on Itokawa will be discussed.

11“Comprehensive analyses of micro-grains from the surface of asteroid returned by Hayabusa:
Implications for collisional space environments”
NAKAMURA, Eizo (Okayama University)
Meteorites are regarded as fragments of asteroids that fall to Earth's surface. However, information
regarding the outer surface of asteroids is presumably destroyed during atmospheric entry, resulting
in our inability to study the solar space-exposed exteriors of planetary bodies other than the moon
sampled by the Apollo missions. JAXA conducted the Hayabusa mission. The original plan of themission was to collect rocks from an S-type asteroid, Itokawa's surface by an impact sampling
method; however, at the time of the touchdown, no projectile was fired, resulting in only minimal
sample recovery from very near the surface. The sample capsule was successfully returned to Earth.
Five lithic grains with diameters near 50 µm, collected during the Hayabusa sample-return mission,
were comprehensively examined for surface texture and geochemistry. Oxygen isotope
compositions indicate that the grains are extraterrestrial. Also considering the major-element
compositions of the grains, it appears that Itokawa's surface is dominantly of ordinary-chondrite
composition. Because the grains were sampled from very near the surface of the asteroid, grainsurfaces retain textures reflecting the space environment influencing the physical nature of theasteroid exterior. The surfaces are dominated by fractures, and the fracture planes contain
sub-um-sized craters and a large number of sub-µm- to several-µm-sized adhered particles, some ofthem glasses. The size distribution and chemical compositions of randomly sampled adhered objects,
and the occurrences of sub-um sized craters, suggest formation by hypervelocity collisions of
micrometeorites down to nm-scale, a process expected in the physically hostile environment at anasteroid's surface.

“Levitation dust: a source of IDPs”
SENSHU, Hiroki (Chiba Institute of Technology Planetary Exploration Research Center)In this talk we will explain the mechanism of photoelectric dust levitation. A dust smaller than 1 µm
can escape the gravitational field of asteroid due to electric repulsion making one of sources of IDPs.

12“An synthesis experiment of GEMS analogue material produced by thermal plasma”
MATSUNO, Junya (Osaka University)
GEMS (Glass with Embedded Metal and Sulfide), which is included in cometary IDPs, is a
sub-micron sized grain composed of an amorphous silicate with embedded Fe, Ni metals andsulfides. Keller and Messenger (2011) proposed that GEMS is condensed products from gases inthe early solar system. However, no condensate experiments reported producing similar materials to
GEMS. ITP (induction thermal plasma) provides so high temperature (~ 10,000 K) environment with
rapid cooling rate (105-106 K/sec) that it makes vapor phase easily and produces nano-sized andnon-equilibrium materials. We injected MgO, Fe, SiO2 on the mean composition of GEMS into ITP.
Run products were characterized by XRD and TEM. Sub-micron sized (about 50 nm) sphericalgrains were observed and their morphology were that each grains have amorphous silicate withembedded an iron (about 10 nm).

“Coagulation and Growth of Dust in Plasma by Charge Fluctuation ~ from Molecules to Clusters
and Dust”
NUNOMURA, Shota (National Institute of Advanced Industrial Science and Technology)
N/A“Coulomb-Crystallization and Behavior of Fine Particles in Experimental Plasmas”
HAYASHI, Yasuaki (Kyoto Institute of Technology Faculty of Engineering and Design)
Because the behavior of fine particles in a experimental plasma on the ground can be directly
observed, the results obtained by observation more or less may contribute to the study on dusts inspace, although conditions of plasmas are different to each other.

Mono-disperse spherical fineparticles can form a Coulomb crystal when they exist in a glow-discharge plasma with a large
amount. We present Coulomb crystals formed by growing fine particles in plasma and therevolutional and rotational behaviors of fine particles in relation to magnetic field.

“Numerical simulation of collisions between dust aggregates consisting of particles with a size
distribution”
WADA, Koji (Chiba Institute of Technology Planetary Exploration Research Center)
We have been performing numerical simulations of collisions between dust aggregates consisting ofequal-sized particles. Particles in proto-planetary disks are, however, expected to have a size
distribution. In this study we carry out numerical simulations of collisions between dust aggregatesconsisting of particles with a size distribution to evaluate the influence of the size distribution on the
dust growth process. A preliminary result will be shown and discussed.

13“Necessary condition for gas giant formation”
KOBAYASHI, Hiroshi (Friedrich Schiller University Jena)
In final stage of planet formation, planetary embryos grow through collisions with planetesimals.
The gravitational scattering by the large embryo induces destructive collisions between
planetesimals. Fragments produced by such collisions are ground down by successive collisions
until the collisional velocities of very small bodies are sufficiently dumped by gas drag. They areconsequently removed by the gas drag.
embryos decreases.

As a result, the surface density of bodies around planetary

Since embryos grow through collisions with such bodies, their growth halts bycollisional fragmentation.

The final planetary mass, which is determined by the equilibrium between

the growth of the embryos and the depletion of planetesimals by collisional fragmentation, is smaller
than the critical core mass for gas giant formation through core accretion. However, since planetaryatmosphere enhance their collisional cross section with small bodies, planetary embryo can exceed
the critical core mass from large initial planetesimals in massive proto-planetary disks.

“Ionization Degrees with Dust Grains in Circum-planetary Disks”
FUJII, Yuri (Nagoya University)
During the formation of gas giants, disks appear around them. They are called circum-planetarydisks. The most promising mechanism of accretion is the turbulence driven by magneto-rotational
instability (MRI). MRI occurs when the ionization degree is high enough for magnetic fields to couple
to gas. We calculate the ionization degrees to investigate if MRI really causes accretion of
circum-planetary disks. Our result suggests that ionization degrees in circum-planetary disks are too
low for MRI to be active.

“Hydrogen ordering in ice observed from neutron diffraction: Application to planetary science”
ARAKAWA, Masashi (Kyushu University)
Laboratory experiments and infrared observations suggest that a lot of water ice exists in a
crystalline phase in our solar system. At ambient pressure, crystalline ice has two kinds of structure,
ice Ih and ice XI. Ice Ih is normal ice with disordered hydrogen atoms. In contrast, ice XI is ahydrogen-ordered phase of ice Ih, and the ordered arrangement of hydrogen atoms makes ice XI
ferroelectric. Existence of hydrogen-ordered ice in space is the subject of continuing astronomical
debate. Electrostatic forces, caused by the ferro-electricity, increase the sticking probability of icy
grains, and might play an important role in grain evolution and planetary formation in space. To
discuss the existence of hydrogen-ordered ice in space, we need to investigate the kinetics andformation process of ice XI. For precise understanding of nucleation and growth process of ice XI, we
performed neutron diffraction measurements of ice doped with various kinds of dopant as a catalyst

14

for the transformation to ice XI. In this presentation, I will talk about how temperature history, kind of
dopant and concentration of dopant affect the formation of ice XI. In particular, I will discuss "memory
effect" of hydrogen ordering in ice.

“Evaporation of icy planetesimals due to planetesimal bow shocks”
TANAKA, Kyoko (Hokkaido University)
We investigated heating and evaporation of a planetesimal due to planetesimal bow shocks
associated with a planetesimal orbiting with supersonic velocity relative to the gas in a protoplanetary
disk. We evaluated surface temperatures and evaporation rates of planetesimals based on the
results of the previous studies of planetesimal bow shocks and formation of protoplanets. We foundthat icy planetesimals evaporate efficiently by bow shock heating in the stage of formation of
protoplanets, where strong bow shocks are produced by gravitational perturbations by the
protoplanets.

“Structure of the Wind-Wind Collision Shocks in Massive Binaries”
OKAZAKI, Atsuo (Hokkai-Gakuen University)
In binaries consisting of two massive stars, strong stellar winds from these stars collide each other.
Shocked region produced by the wind-wind collision is thought to be a place where not only theparticle acceleration and resulting high-energy emission but also the dust formation takes place.
Studying the wind-wind collision shocks is thus of multi-facet importance. Numerical simulations suit
it very well, given that massive binaries, in general, have large orbital eccentricities. In this talk, I willreport on the structure of the wind-wind collision shocks in two massive binaries, eta Carinae andWR 140, obtained from 3-D numerical simulations.

“Interactions between stellar winds and the interstellar medium and their effects on dust grains”
UETA, Toshiya (University of Denver)Latest results of far-IR imaging by AKARI, Spitzer, and Herschel have suggested that interactions
between stellar winds and the interstellar medium (ISM) are not uncommon. These findings imply
that circumstellar dust grains can be processed at the interface between stellar winds and the ISM
before they are injected into the ISM and become ISM dust grains, and that there may be significant
implications to our understanding of "Galactic dust grains". I will review the latest far-IR observations
on the interface regions between stellar winds and the ISM and discuss what follow-up studies with
ALMA and SPICA will help enhance our understanding of dust "production" by stars.

15“Properties of dust and PAHs in clouds with high-energetic phenomena based on the AKARI
all-sky-survey maps”
KONDO, Toru (Nagoya University)
Very recent radio, X-ray, and gamma-ray observations have revealed various kinds of high energetic
phenomena from interstellar clouds in our Galaxy, such as the shock heating of clouds bysupernovae and the production of TeV gamma-rays through interaction of high-energy cosmic-ray
particles with clouds. Under such harsh interstellar environments, PAHs are processed anddestroyed quite easily, while submicron dust grains survive for a relatively long time, which arecollisionally heated to high temperatures. Therefore, spatial variations of the properties of dust andPAHs can provide observational evidence that energetic phenomena really take place in thecorresponding regions from entirely different aspects from the radio, X-ray, and gamma-ray
observations. We are currently conducting intensive reduction of AKARI all-sky-survey diffuse data,
aiming at the public release of the all-sky maps within a year. By using part of the diffuse maps so far
created at 9, 18, 90, and 140 micron, we have made a correlation study of the AKARI all-sky survey
with the NANTEN, Suzaku, and HESS data. Among the AKARI data, the 9 µm map is sensitive only
to PAHs without contribution from the mid-IR continuum, while the 18 micron map is important totrace shock-heated hot dust. As a result, we find that the ratios of the PAH to dust emission
brightness significantly change in the clouds likely associated with high energetic phenomena. In thepresentation, we summarize the overall results so far obtained, and also briefly introduce ouractivities on the all-sky-survey data reduction.

“Distribution of the 3.3µm PAH band in our Galaxy detected by the AKARI IRC”
TSUMURA, Kohji (ISAS/JAXA)
Polycyclic aromatic hydrocarbon (PAH) is a good tracer of distribution and physical properties of theinterstellar dust because PAH has distinct emission bands in the near-infrared region. We analyzed
the spectral data of the interstellar diffuse emission obtained with the low-resolution prism
spectroscopy mode on the AKARI Infra-Red Camera (IRC) NIR channel in (2.0-5.0µm) wavelength
region.

As a result, a good correlation of the 3.3µm PAH band emission to the thermal emission

from the interstellar dust was confirmed. We also construct the reference diffuse Galactic light (DGL)
spectrum at 2-5µm by using this 3.3µm PAH band correlation as a tracer.

16

Proceedings

17

A study of dust grain formation: A role of minor elements
Akemi Tamanai1,2, Harald Mutschke2, Jürgen Blum3
1. Kirchhoff Institute for Physics, Heidelberg University, Im Neuenheimer Feld 227, D-69120 Heidelgerg,
Germany
2. Astrophysical Institute and University Observatory, Friedrich-Schiller-University Jena, Schillergaesschen 3,
D-07745, Germany
3. Institute for Geophysics and Extraterrestrial Physics, Technical University Braunschweig, Mendelssohnstr. 3,
D-38106, Braunschweig, Germany

Abstract
A large amount of silicate dust grains has been detected in many astronomical environments.
However, silicate dust grain formation itself is not well understood yet. Although theoreticalsimulations for possible condensates illustrate how abundant the different silicate dust materials
are in astronomical surroundings and providing us with a key to clarify the formation of silicates,
these simulations do not cover all dust species which possibly exist. In order to deepen theunderstanding from the first species condensing out of the gas to the entire condensation
sequence of dust materials present in the environments, we take into account the species whichcontain minor elements (Al, Na, and Li) and perform the mid-IR spectroscopic measurements for
obtaining benchmark spectra.
1. Introduction
Olivine as well as pyroxene dust grains have been observationally confirmed in various
astronomical environments. Still it is not straightforward to obtain the precise chemicalcompositions (e.g. the ratio of Mg and Fe) and morphological aspects (e.g. size, shape,
aggregate state) through the observed spectra. Meteorites, on the other hand, have provided us
with exactly these detailed information. For instance, more than 10 different types of silicateshave been discovered in the Allende meteorite including phyllosilicates (e.g. Mg 3Si4O10(OH)2(talc), (Mg, Fe)7Si8O22(OH)2 (anthophyllite)), but also different types of oxides and fullerene-likecarbon (Jarosewich et al. 1987; Brearley 1996; Rubin 1997; Harris et al. 2000; Nozawa et al.
2009; Ma et al. 2010). Although it is possible to extract and analyze a small amount of a particular
dust grain embedded in a meteorite, it is difficult to identify that via observations if the dust
species is small in quantity. In other words, only very abundant dust grains can be detected andshow their absorption peaks in the mid-infrared (MIR) regions of the electromagnetic spectrum.
Though the less abundant grains definitely take part in the dust evolution process, the chemical
pathways with regard to dust grain formation have not clearly understood yet. Especially, theformation process of one of the most abundant dust grains, i.e. the silicates, has not been
elucidated evidently so far so that questions like the following remain: What is the highesttemperature condensation material? How does silicate dust formation take place? How muchsilicates are formed? Does silicate formation start by surface growth on high-temperaturecondensates such as, corundum and spinel, or directly from the gas phase? Understanding theevolution and formation of silicate dust grains can lead us not only to the properties andprocessing history of the entire grain population, but also to the transition and role of them invarious surroundings in connection with star and planet formation. For that, we have taken into
account the less abundant elements (Na, Li, Al), compared to Si and Mg, and have performed
MIR spectroscopic extinction measurements of Al2SiO5 (kyanite), NaAlSi2O6 (jadeite), and
LiAlSi2O6 (spodumene) by means of aerosol spectroscopy. The outcoming spectra will set criteriafor MIR spectra obtained by condensation experiments in our group in the future (see moredetails in Wetzel et al. in this proceedings).

18

Polyester membrane filter

FTIR spectrometer
Bruker 113 v
Muti-path gas cell

DTGS Detector

Incident light
Parabolic mirror

Flat mirror

Parabolic mirror

MirrorMirrorImpactor
Pump

Dust flow
generator

Fig. 1 A schematic diagram ofthe experimental apparatus foraerosol spectroscopy.

2. Experimental Setup
Fig. 1 shows the setup of the apparatus for the aerosol spectroscopic measurement. We have
utilized a dust flow generator (Palas RBG1000) for dispersing a powdered sample into a nitrogen(N2) gas stream. The dense aerosol is carried to a two-stage impactor which separated the large
particles from the small ones (< 1 m). The small-sized grains are concentrated by the impactor
until a concentration of 106 particles per cubic centimeter is reached. In the end, only the small
particles arrive in a White-type long-path infrared cell (MARS-8L/20V, Gemini Scientific Instr.)
which retained 18 m path length by taking advantage of multiple reflection between two gold
mirrors mounted on both sides in order to intensify the sensitivity. Since a Fourier Transformation
Infrared Spectrometer (FTIR, Bruker 113v) with DTGS detector with CsI windows is fixed to thecell, it is possible to measure the extinction spectrum of the suspended dust particles in N 2 gas.
As a consequence, we can avoid the environmental effect during the aerosol measurements, andthe measurement condition is closer to a vacuum with regard to the dielectric function of medium(N2: 1.0) compared with the KBr (potassium bromide; =2.3). Moreover, the risk of causing thepowdered sample structurally deforms by the high pressurization required for the pellet technique,this problem is avoidable with aerosol spectroscopy. For clarifying the environmental effect, wehave performed extinction measurements obtained by CsI (cesium iodine) pellets as well. CsI isparticularly transparent between 2 and 50 m in wavelengths, and has a higher value of opticalconstants than that of KBr (CsI n=1.74 & KBr n=1.52). The aerosol particles are extracted on apolyester-membrane filter that is mounted a little above the cell. The morphological aspect of theextracted particles was investigated by making use of a scanning electron microscope (SEM)(see more details in Tamanai & Mutschke 2010).

3. Samples
Silicate minerals, which contain the less abundant refractory elements, such as Al, Na, and Li, arediscussed in this investigation.
A chemical composition of Al2SiO5 is a polymorph with three minerals, which are kyanite,
andalusite, and sillimanite. Kyanite is formed under conditions of low temperatures (< 800°C) andhigh pressures, whereas sillimanite formation takes place at high temperatures (> 800°C) and a
relatively wide range of pressures. On the other hand, andalusite is formed at a limited range oftemperatures (~200°C and < 800°C) and pressures (< 4 kbar) (e.g. Althaus 1967; Whitney 2002).
Here we only consider kyanite, which belongs to triclinic crystal system (e.g. Deer et al. 1997).

19

The presence of corundum (Al2O3) and quartz (SiO2) induce formation of kyanite. This reversible
chemical reaction will be described as(Mueller et al. 1977).
According to the equilibrium condensation calculations for the solar nebular (Davis & Richter2003), albite (NaAlSi3O8) formation is expected at temperatures of approximately < 1000K andpressure of 10-3 bar. In fact, albite has been found in a unique unequilibrated chondrite,
Kakangari (Brearley 1989) and a carbonaceous chondrite (type CR2), the Kaidun meteorite(Ivanov et al. 2001). Albite is not a target material for this investigation. It is closely related tojadeite (NaAlSi2O6) which has a monoclinic crystal structure and is possibly formed by a chemical
reaction of jadeite and quartz. Since this is a reversible chemical reaction as well, albite
decomposes to yield jadeite and quartz as below(Bucher & Grapes 2011).
Formation of spodumene (LiAlSi2O6: lithium aluminum silicate) occurs at low temperatures(<700°) and high pressure (> 2 kbar) (Stewart 1978) and has a monoclinic crystal structure (e.g.
Deer et al. 2004). These rare elements have been engaged in discussions of isotopes inmeteorites. Hydrogen, Helium, and Lithium were produced in the big bang nucleosynthesis sothat Lithium is one of the primordial elements. According to Israelian et al. 2009, Lithium has been
depleted effectively in our Sun like planet-hosting stars as compared with planet-free ones. In thismanner, the existence of Lithium plays an important role in extrasolar planet search as well as
mineral formation with different lithium isotopes.
Unfortunately, kyanite, jadeite, and spodumene have not been detected observationally, nor
in meteorites. However, it is still worth investigating the optical properties of these materials for
future potential. Table 1 represents a list of the measured samples, including their chemical andphysical properties. All three samples are naturally formed minerals. Thus, the relative
abundance of the elements (subscripts) is not exactly like the chemical formulae in Table 1. Notethat fact-based chemical formulae (by EDX analysis) will be described elsewhere.

MineralogicalName
KyaniteJadeiteSpodumene

Table 1: Properties of the three samples
Chemical Formula
Product
Particle size(m)
Al2SiO5Natural<1
NaAlSi2O6
Natural (China)<1
LiAlSi2O6
Natural (South Africa)<1

Particle shape(m)
Irregular
Irregular
Irregular

4. Extinction Measurements
As described in Sec. 2, we have performed extinction measurements of kyanite, jadeite, andspodumene by making use of aerosol spectroscopy in the wavelength range between 2 and 50
m. In order to clarify the environmental effect by the surrounding medium, spectra obtained by
the aerosol spectroscopy are plotted together with spectra taken by the CsI pellet technique in Fig.
2. All plotted spectra are normalized since aerosol spectroscopic measurements are notquantitative. Redshift of bands can be observed along the entire extinction spectra of the threematerials investigated by the pellet technique, due to the effect of electromagnetic polarizationinduced by the medium, CsI (e.g. Tamanai et al. 2006). Also, the disparity in band profiles

20

between the aerosol and CsI pellet measurements may result from: (1) sample structure
deformation due to the high pressurization for making the CsI pellets, (2) the homogeneity anduniformity of particle dispersion in the measurement environment (localized or delocalized), (3)particle orientation, (4) distortion of crystal lattice during the grinding procedure (by hand or by aball mill) (Imai et al. 2009), (5) the particle morphology, especially size and shape of aggregates(see also Tamanai et al. 2009).
1.2(a) Kyanite

CsI
AerosolNormalized Extinction

1
0.8
0.6
0.4
0.2
0
8

9 10

20
Log wavelength [m]

30

40

Aero.
9.88
16.31
22.55
28.32
36.14
CsI
9.69
14.85
18.32
25.96
32.28
42.40

10.49
17.60
23.90
29.22
39.03

11.10
18.27
24.88
30.39
40.71

13.66
19.27
25.25
31.40
42.36

14.64
20.34
25.91
32.22

15.64
21.23
26.63
33.79

10.03
15.70
19.81
26.80
33.82

10.59
16.06
21.42
28.36
36.08

11.16
16.52
22.88
29.25
37.72

12.17
17.08
24.91
30.44
39.09

14.01
17.60
25.36
31.46
40.77

9.34
15.04
26.74

9.52
16.71
30.00

9.86
19.25
31.75

10.61
21.13
34.42

11.65
22.77
38.75

9.42
15.07
27.24

9.52
17.04
30.17

10.02
20.22
31.85

10.86
21.60
34.47

11.67
22.84
38.96

9.16
19.69
32.66

9.56
20.65
35.65

10.69
24.85
42.32

11.60
26.33

15.44
29.38

9.22
18.76
36.34

9.72
20.53
42.69

10.95
26.70

11.61
29.68

15.66
31.68

1.2
(b) Jadeite

CsI
AerosolNormalized Extinction

1

Aero.
8.94
13.37
25.01
CsI
9.19
13.39
25.11

0.8
0.6
0.4
0.2
0
8

9 10

20
Log wavelength [m]

30

40

1.2
(c) Spodumene

CsI
AerosolNormalized Extinction

1
0.8
0.6
0.4
0.2

Aero.
8.65
16.59
31.25
CsI
8.93
16.99
32.79

0
8

9 10

20
Log wavelength [m]30

40

Fig. 2 Normalized extinction vs. log wavelength of (a) Kyanite; (b) Jadeite; (c) Spodumene
spectra obtained by the aerosol and CsI pellet measurements. The black solid line indicates theaerosol measurements, and the grey dotted line denotes the CsI pellet measurements. All peak
positions are listed in the tables on the right hand side. The boldface numbers are the strong
peaks in each material.

21

Both Jadeite and spodumene belong to inosilicates, which can be expressed by a general
structural formula of XYZ2O6. In the case of jadeite and spodumene, the octahedral site (X) is
filled with a different element (Li or Na), whereas Y and Z are occupied by the same elements(Y=Al and Y=Si). However, the strongest Si-O stretching vibration band of jadeite is situated to
the long wavelength side (9.56 m) compared to that of spodumene (9.34 m).
5. Summary and Outlook
MIR spectroscopic measurements of silicate materials (kyanite, jadeite, spodumene), whichcontain less abundant elements such as Al, Na, and Li as compared with Si and Mg, areperformed in medium-effect-free environments. These experimentally obtained spectra will bebenchmarks for spectr, which will be measured by condensation experiments in the future.
Acknowledgements
We greatly appreciate Prof. A. Pucci for giving us an opportunity to present this project at GFW. We
are thankful to Prof. H.-P. Gail for insightful opinions and comments. We are grateful to G. Born for her
assistance with sample preparation and analysis and W. Teuschel for his technical support of ourexperimental apparatus. A. Tamanai expresses to Prof. T. Yamamoto her gratitude for the financial
support of the workshop participation and CPS visit. Astronomical research has been supported by
Deutsche Forschungsgemeinschaft (DFG) under the grant FOR 759.
Reference:
Althaus, E. 1967, Contr.Mineral.and Petrol., 16, 29-44
Brearley, A.J. 1989, Geochim. Cosmochim. Ac., 53, 2395-2411
Davis, A.M., and Richter, F.M. 2003, in Meteorites, Comets, and Planets, ed. A.J. Davis, Vol.1,
Treatise on Geochemistry, eds. H.D. Holland, and K.K. Turekian (Oxford: Elsevier-Pergamon)
Brearley, A.J., and Prinz, M. 1996, LPS, 1081, 161-162
Bucher, K., and Grapes, R. 2011, Petrogenesis of Metamorphic Rocks (Berlin: Springer-Verlag)
Deer, W.A., Howie, R.A., and Zussman, J. 1997, Rock-Forming Minerals: Orthosilicates (UK: The
Geological Society)
Deer, W.A., Howie, R.A., and Wise, V.S. 2004, Rock-Forming Minerals: Framework Silicates, Silica
Minerals, Feldspathoids, and Zeolites (UK: The Geological Society)
Harris, P.J.F., Vis, R.D., and Heymann, D. 2000, Earth and Planetary Science Letters, 183, 355-359
Imai, Y., Koike, C., Chihara, H., Murata, K., Aoki, T., Tsuchiyama, A. 2009, A&A, 507, 277-281
Israelian, G., Mena, E.D., Santos, N.C., Sousa, S.G., Mayor, M. et al. 2009, Nature, 462, 189-191
Ivanov, A.V., Kononkova, N.N., Zolensky, M.E., Migdisova, L.F., and Stronganov, I.A. 2001, Lunar andPlanetary Science, XXXII, 1386
Jarosewich, E., Clarke, R.S.Jr., and Barrows, J.N. 1987, The Allende Meteorite Reference Sample(Washington D.C.: Smithsonian Institution Press)
Ma, C. 2010, American Mineralogist L., 95, 188-191
Mueller, R.F., and Saxena, S.K. 1977, Chemical Petrology with applications to the Terrestrial Planets
and Meteorites (N.Y.: Springer-Verlag)
Nozawa, J., Tsukamoto, K., Kobatake, H., Yamada, J., Satoh, H. et al. 2009, Icarus, 204, 681-686
Rubin, A.E. 1997, Meteoritics & Planetary Science, 32, 231-247
Stewart, D.B. 1978, Am. Mineral., 63, 970-980
Tamanai, A., Mutschke, H., Blum, J., and Meeus, G. 2006, ApJ, 648, L147-L150
Tamanai, A., Mutschke, H., Blum, J., Posch, Th., Koike, C., and Ferguson, J.W. 2009, A&A, 501, 251267
Tamanai, A., and Mutschke, H. 2010, Infrared spectroscopy of dust particles in aerosols for
astronomical application, in Fundamentals and Applications in Aerosol Spectroscopy, eds. J.
Reid and R. Signorell, pp.101-123 (Florida: Taylor & Francis Group)
Whitney, D.L. 2002, Am. Mineral., 87, 405-416

22

IR spectra of silica (SiO2 ) polymorphs
Koike C., Imai Y., Noguchi R., Chihara H.,
Ohtaka O., Tsuchiyama A.,
Department of Earth and Space Science,
Graduate School of Science, Osaka University,

Abstract
Recently, silica SiO2 was detected in the protoplanetary disk around the T Tauri stars and in
the debris disks. In some objects, annealed silica (cristobalite and/or tridymite) was main mineralamong observation. Contrary to the silica, Mg-rich silicates such as forsterite (Mg2 SiO4 ) and enstatite(MgSiO3 ) were detected as for main mineral of circumstellar dust. This referred silica is quartz or
annealed silica among the various crystal type such as quartz, fused quartz, cristobalite, coesite,stishovite. We prepared these various type of silica in natural, synthetic, and commercial. We reporton the infrared spectra of these silica and discuss on the possibility of detection silica for observation.

1

Introduction

Silica (silicon dioxide SiO 2 ) of terminology is used in solid state and astronomy as for fused quartz,
vitreous silica, quartz, cristobalite, etc. in case by case, and many workers were confused to use thisterminology in many ﬁeld. Especially, astronomers used Silica of terminology as for quartz in some case,and for annealed silica in some case. Astronomers used Silica not giving a clear deﬁnition of this term.
Mineral of Silica was taken into as major component to explain the spectra of T Tauri star HEN
3-600A, and in this case Silica is introduced into as quartz (Honda et al. 2003). A few T Tauri stars(TTS) show prominent narrow emission features indicating annealed silica, where two polymorphs ofsilica, tridymite and cristobalite, are the dominant forms in the TTS (Sargent et al., 2009). The annealed
silica was ﬁrstly introduced into the component of circumstellar dust, and this ’annealed silica’ wassynthesized after annealed of fused quartz at 1220 K (947 ◦ C) for 4.5-5 h (Fabian et al. 2000). Recently,
some strange features were detected in debris disks (HD 172555, HD 15407A, HD 23514A, HD 23514),
where prominent peaks showed at about 9.3, 12, 21 μm (Lisse et al. 2009), and at about 9, 16, and20-22 μm (Melis et al. 2010, Rhee et al., 2008, Fujiwara et al., 2009). HD 172555 has infrared spectra
at about 9.3 and 19 μm with subtle peaks at 7.5 - 8.2, 11.2 , and 33 μm, where the spectra of silica suchas obsidian and tektite were well ﬁtted to the observation (Lisse et al., 2009). The strange spectra of
HD 15407 by Fujiwara et al.(2009 (GFWS)) were noticed especially and these spectra were ﬁtted withamorphous silicate (amorphous pyroxene) and fused silica. But, the 16 μm band was not ﬁtted by theirmodel (Fujiwara et al. 2009).
Further, coesite and stishovite were discovered in lunar meteorite Asuka-881757 (Ohtani et al. 2010),
and these high-pressure polymorphs of SiO2 indicate that the meteorite experienced an equilibrium
shock-pressure of at least 8-30 GPa. High-pressure impact metamorphism and formation of high-pressureminerals may be common phenomena in debris disks.
These minerals of silica polymorphs, if exist in circumstellar or debris disks, may be ascertainedby those of infrared spectra. Neverthless, the infrared spectra of silica polymorphs in wide wavelength
region were not measured in laboratory, where even measured data of silica polymorphs were only limited
wavelength region.
Here we measured the infrared spectra of the various type of silica in wide wavelength region. Following
to our result for annealing fused quartz, we discuss on the possibility of detection silica for observation.
Proceedings of 29th Grain Formation Workshop 2011, Nov.9

1

23

Figure 1: Infrared spectra of silica for present samples.

2

Preparing samples and measurements

We prepared many kinds of silica such as fused quartz, quartz, cristobalite, coesite and stishovite.
• Fused quartz
Commercial fused quartz were obtained by Nippon Aerosil Co. Ltd. (details in Koike et al. 1994).
• Quartz
Natural quartz (production area is not known) was ground by hand and further ground by ball mill
gently.
• Cristobalite
Fine powder of cristobalite were prepared by synthesized, by commercial, and by natural samples.
Fine powder of fused quartz were annealed at 800 - 1200 ◦ C for various duration times, and 1500 ◦ C
for 3 and 6 days. Further, ﬁne powder of commercial cristobalite was obtained from Omura Ceratec
Inc.(Japan). Natural cristobalite was obtained from Cougar, Siskiyou Co., California, USA.
• Coesite
The ﬁne powder of fused quartz was pressed at 6 GPa and 1300 ◦ C for 30 min at Osaka University.
The compressed powder became coesite, which was conﬁrmed with XRD.
• Stishovite
The natural quarzt was compressed at 13 GPa and 1200 ◦ C. The samples was synthesized by Tani
(2000), and kindly given to us for measuring infrared spectra.
2

24

Figure 2: infrared spectra of cristobalite and coesite compared with that of forsterite.
The ﬁne particles of silica polymorphs were well dispersed in KBr and polyethylen powder, and pressed
into KBr pellts and polyethylene sheet. Then the transmittance of pellets and sheets were measured with
Nicolet 670 and/or 6700 in mid- and far-infrared region. The mass coeﬃciency κ was derived fromtransmittance T as follows.
S
1
κ=
ln
(1)
M T
Here, S is surface area of pellets or sheets and M is weight mass of sample particles in pellets or sheets.

3

Results and discussion

3.1

Cristobalite

The powder of fused quartz was annealed at 1000 ◦ C for 1 h following to the results by Fabian et al.,
where cristobalite appeared at 1220 K (947 ◦ C) for 4.5 - 5h. But, these samples had no change of infrared
spectra from starting spectra. Further changing the temperature and duration time for annealing fusedquartz at 800-1000 ◦ C for 1 - 7h, the infrared spectra were also not changed from the starting spectra.
Even after annealing at 1100, 1150 , 1200 ◦ C and for long duration time (2 days - 8 days), the fusedquartz was not changed into cristobalite with IR spectra, although small quantity of cristobalite wasascertained by XRD. This is inconsistent to the results by Fabian et al. (2000), where fused silica wasannealed at 1220 K (947 ◦ C) for 4.5-5 h and changed into cristobalite both by IR and XRD. Usually insolid state ﬁeld, cristobalite appeared at above 1400 ◦ C by phase diagram of silica ( Šimon et al. 1953,
Plendl et al. 1967, Swaison et al.2003, Lakshtanov et al. 2007). Following to phase diagram of silica, the3

25

Figure 3: The infrared spectra of coesite and stishovite (?) . The sample of stishovite is not conﬁrmed
and marked as ?.
cristobalite changed after annealing at low temperature of about 950 ◦ C is very curious result, and thismay be due to some impurity eﬀects or some unique condition.
After annealed fused quartz at 1500 ◦ C for 3 and 6 days, cristobalite was conﬁrmed with XRD, andthe infrared spectra were measured to be same with each other, and similar with the spectra by Fabian
et al.(2000). Where some peaks were shifted, and the peak of 68 μm band was shifted about 1 μm. The
infrared spectra of annealed silica at 1500 ◦ C were almost same with those of commercial cristobalite.
Natural cristobalite showed strange spectra due to contamination of some minerals.
Here after, the infrared spectra of cristobalite were shown as annealed at 1500 ◦ C for 3d. The
characteristic peak appeared as the strong and broad band at 68 μm. The 68 μm band show very broadband compared with the 69 μm band of forsterite (Mg2 SiO4 ). But, the 33 μm band of cristobalite showedvery sharp compared with the 33 μm band of forsterite. Forsterite is main component of circumstellar
dust and observed as prominent peaks at 11.2, 16, 33, and 69 μm.

3.2

Coesite

The infrared spectra of coesite were shown as many sharp ﬁne peaks in Fig. 1 and 2, and as the sharp
peak detected at 69 μm, which is very similar to those of forsterite. The 69 μm band of coesite is very
sharp compared with that of forsterite. This peak was ﬁrstly detected deﬁnitely.

3.3

Stishovite

The spectra of stishovite showed similar to those of coesite. Usually, stishovite showed strong peaks at
about 12-18 μm compared with the data by Lyon (1962), and the present spectra is very curious. So, we
marked this sample as ? until it is conﬁrmed to ascertain as stishovite. The similar spectra as coesite
might be due to some aging eﬀects, but it is not clear. We must make clear why this strange spectra
were appeared.

3.4

Fused quartz and Quartz

The spectra of fused quartz (400 Åsample; average size is about 40 nm) were shown in Fig. 1 (Koike etal. 1994). The three broad peaks appeared at 8.9, 12.3, and 20.8 μm.
The spectra of quartz showed sharp peaks at about 9, 12, 14, 19-21, 25, 27, and 38 μm. The weak
peak at 78 μm was very subtle, and was ascertained with bulk quartz sample (Russel et al. 1967). This
peak was not appeared with reﬂectance measurements due to weak band by FT/IR Nocolet 670 , and4

26

the measured reﬂectance spectra between 270 - 2000 cm−1 were almost similar as those by Spitzer andKleinman (1961).

3.5

All samples

The spectra of all samples are compared with each other. The strength and peak position of 9 μm band
are almost all similar for all samples. But, for another bands, many sharp peaks appear for coesite andquartz. In far infrared region at about 69 μm, the strong and broad band appears for cristobalite, andthe strong and sharp peak appears for coesite.

3.6

Spectra at low temperature compared with those of forsterite

The spectra of cristobalite and coesite show similar bands at 16, 33 and 69 μm as those of forsterite. The
mineral of forsterite is main component in circumstellar dust, and ascertained by peak positions from
infrared spectra of observation. In order to discriminate between the samples (cristobalite and coesite)
and forsterite, the spectra of coesite and cristobalite were measured at low temperature from 8 K to
RT. Those spectra at low temperature were shown in Figs. 4, and the peak positions of cristobalite andcoesite were compared with those of forsterite as shown in Fig. 5 and 6. Most peak positions of eachband shifted to shorter wavelength as cooling from RT to 8 K, but some peak positions shifted to longer
wavelength as cooling to 8 K. In the case of cristobalite, the peak positions of 33 and 36 μm bands shift
to longer wavelength as cooling to low temperature. In the case of coesite, those of 34, 38, 39, and 69
μm bands shift to longer wavelength as cooling to low temperature. The diﬀerences between cristobalite,
coesite and forsterite are clear for 16 μm band as peak positions are diﬀerent and distinguished between
each other. As for 33 μm band, the peaks of cristobalite and coesite are sharper than that of forsterite.
As for 69 μm band, the band of cristobalite is very broad and peak position is clearly diﬀerent from thatof forsterite. The band of coesite is very sharper than that of forsterite, and shift to longer wavelengththan that of forsterite, although the peak position became same as that of forsterite at 200 K .

Figure 4: The spectra of cristobalite and coesite at low temperature.

5

27

Figure 5: The peak positions of cristobalite and coesite at low temperature compared with those of
forsterite.

Figure 6: The peak position of 69 μm band of cristobalite and coesite at low temperature compared withthose of forsterite.

6

4

Implication to the observation

28

Cristobalite was ascertained with XRD and infrared spectra after annealing of fused quartz at 1500 ◦ C
for 3 days and 6 days. This annealing temperature is consistent with the results by many workers (Šimon
et al. 1953, Plendl et al. 1967, Swaison et al. 2003). The temperature of 1220 K (947 ◦ C) by Fabian
et al. is too low compare with our value and may be due to some unique condition. Around T Tauristars, annealed silica (cristobalite and/or tridymite) would exist from observation of mid infared , whereannealing amorphous silica around T Tauri stars was proposed (Sargent et al. 2009). In this case, it isvery important whether temperature of the circumstellar disk can be raised to above 1400 - 1500 ◦ C or
not instead 947 ◦ C. If annealed temperature could not be raised to above 1400 - 1500 ◦ C, another process
such as collision must be introduced.
As for the evidence of existing cristobalite, the broad and strong band at about 68-62 μm could be
detected if more detailed observation would be worked in far infrared region. If dust mineral of coesite
and cristobalite might exist around debris disks production by high speed collisions among planetesimals,
the charcteristic bands could be detected by the high resolution observations in 60-70 μm region.
As for the observed spectra of the debris disks of HD 15407 (Fujiwara et al. 2009), the present
cristobalite may be ﬁtted well the observed peaks at 9, 16, and 20-22 μm band instead of fused quartz.
The broadness of peaks may be explained after heating up to RT to 250 ◦ C, where 16 μm band exist
clearly (Zhang et al. 2007), instead of increase of size of dust particles.

Acknowledgements
We wish to thank Dr. A. Tani for providing the stishovite samples.

References
Fabian D., Jäger C., Henning T., et al., 2000, A&A, 364, 282
Fujiwara et al., Grain Formation Workshop, 2009, pdf ﬁleHonda M., Kataza H., Okamoto Y., et al. , 2003, APJ, 585, L59
Koike C. and Shibai H., 1994, Mon. Not. R. Astron. Soc., 269, 1011
Lakshtanov D. L. and Sinogeikin S. V., 2007, Phys. Chem. Minerals 34, 11
Lisse C. M., Chen C. H., Wyat M. C., et al., 2009, APJ, 701, 2019
Lyon R. J. P., 1962, Nature 196, 267
Melis C., Zuckerman B., Rhee J., and Song I., 2010, APJ 717, L57
Ohtani E., Ozawa S., Miyahara M. et al. 2010, PNAS 108, 463@
Plendl J. N., Mansur L. C., Hadni A., et al., 1967, J. Phys. Chem. Minerals, 14, 527
Rhee J., Song I., and Zuckerman, 2008, APJ 675, 7777
Russell E. E., and Bell E. E., 1967, J. Optical Society America 57, 341
Sargent B., Forrest W., Tayrien C., et al., 2009, APJ 690, 1193
Šimon I. and McMahon H. O., 1953, J. Chem. Phys. 21, 23
Spitzer W. G. and Kleinman D. A., 1961, Phys. Rev. 121, 1324
Swaison I. P., Dove M. T., and Palmer D. C., 2003, Phys. Chem., Minerals. 30, 353
Tani A., Yamanaka C., Ikeya M., et al., 2000, Appl. Magn. Reson. 18, 559
Zhang M., and Scott J. F., 2007, J. Phys.: Condens. Matter 19, 275201

7

29

Vibrational spectroscopy of SiO and Mg silicate systems
S. Wetzel, A. Tamanai, and A. Pucci
Kirchhoff-Institut für Physik der Universität Heidelberg;
Im Neuenheimer Feld 227; 69120 Heidelberg

Abstract
Silicates are the most abundant minerals in space and their IR spectral features depend on physical properties such as the chemical composition and on theformation conditions e.g. temperature and pressure. Therefore they can provide
information not only on the relative abundance of the elements, but also on formation processes and the environment in which the dust is formed. To perform
sophisticated analysis of IR spectra obatined by space observations a comprehensive database is required, including materials with various compositions andinterfacial interaction effects have to be considered. In order to deepen the understanding of IR spectra of silicates, we performed in situ IR spectroscopy oflayers produced by condensation of evaporated SiO, Mg and Fe, respectively.
Details on the experimental setup and first results from the measurements arepresented in this work.

I.

INTRODCUTION

The strong absorption feature at around 10 µm is attributed to the Si-O stretching vibration which is sensitive to changes of surrounding bonding parameters as it has been for
example demonstrated for the SiO/Si(111) interface1 . To better understand the IR spectraof silicates, we investigate materials involved in silicate formation by in situ IR spectroscopy
under ultra high vacuum (UHV) conditions.
1

30

II.

EXPERIMENTAL

A schematic drawing of the experimental setup is shown in Fig. 1. It consists of a FTIRspectrometer that is coupled to an UHV chamber, with a base pressure of 1 × 10−10 mbar,
in which the condensation experiments were performed. Equipped with a triple evaporator that allows individual, subsequent or simultaneous evaporation of up to three different
materials one can measure single layers, multilayer structures and layers produced by coevaporation. The evaporation rate is determined with a quartz crystal microbalance (QCM),
that can be inserted into the particle beam, before and after each spectroscopic experiment.
With a known density of the material, this allows to calculate the film thickness corresponding to each measured spectrum. The sample holder is attached to a cryostat andcan be heated by electron bombardment, that allows to vary the substrate temperaturefrom 40 to 1300 K; therefore, it also makes annealing experiments possible. The chamber isadditionally equipped with a low-energy-electron-diffraction (LEED) system to investigate
crystalline surfaces, a gas system and quadrupole mass spectrometer (QMS) for experiments
under well-defined partial pressures of gases. The beam path outside the UHV chamber isevacuated to 4 mbar. IR spectra were measured with a mercury-cadmium-telluride (MCT)
detector in the wavelength range between 2 µm (5000 cm−1 ) and 22 µm (450 cm−1 ) with aspectral resolution of 4 cm−1 under normal incidence of light and normalized to the spectrumof the bare substrate. Si(111) wafers with dimensions of 10 × 10 mm2 and a thickness of
0.525 mm were used as substrate. Before deposition of the metal films the natural oxide wasremoved by heating the wafer to 1300 K. In the case of SiO deposition the natural oxide wasleft on the surface as it only effects the first stages of film growth below 1 nm film thickness.1

III.

BASICS

Approximate transmittance formulas of a thin film with the thickness d (d  λ (wavelength)) can be derived by using the Fresnel equations.2 On the other hand, the normaltransmittance of a thin film on a thick substrate (d  λ)3,4 is given by
Trel ≈ 1 −2d ωc√ Im(f )1 + s

2

(1)

31

UHV chamber(base pressure < 10-10 mbar)
gas system

LEEDFTIR-spectrometer(Bruker IFS66v)

sample attachedto cryostat/heater(40 - 1300 K)

IR beamdetector(MCT/DTGS)

triple evaporator

quartz crystal
microbalance (QCM)
quadrupole mass
spectrometer (QMS)

FIG. 1: A schematic drawing of the experimental setup. A FTIR spectrometer (Bruker IFS66/v)
is coupled to an UHV chamber with a base pressure of 1 × 10−10 mbar. The setup is equipped witha triple evaporator that allows individual, subsequent or simultaneous evaporation of up to threedifferent materials. The temperature of the substrate can be controlled between 40 and 1300 K(ω = 2π λc ; s : the dielectric function of the substrate; f : the dielectric function of thefilm). Only the transverse optical (TO) modes with their frequency at the maximum of theimaginary part of the dielectric function  can be detected under normal incidence of light .
An appropriate model has to be used to determine the dielectric function from such ameasurement.5 The oscillator model introduced by Brendel6
NX

1
(ω) = ∞ +
2πσj
j=1

Z ∞

−

e
−∞

(x−ω0,j )2
2σ 2
j

Ωj
dx ,
2
x − ω 2 − iωγj(2)

has been used to describe our measured spectra, as it is especially suited to describe disordered amorphous systems. It consists of a dielectric background ∞ and N Gaussian
distributions of Lorentz oscillators with resonance frequencies ω0,j , damping constants γj ,standard deviations σj of the Gaussian distributions and oscillator strengths Ωj . For are
more detailed description we refer to the work by Klevenz et. al.5 and the references therein.

IV.

RESULTS & DISCUSSION

Relative transmittance spectra taken during the condensation of SiO on a Si(111) substrate at room temperature are shown in Fig. 3 . Two broad peaks with resonance frequencies
at 10.163 µm as well as 14.085 µm and a shoulder at longer wavelength develop as film thick3

32

1 4 . 0 8 5 µm

1 0 . 1 6 3 µm

0 n m

re l. tra n s m itta n c e

1 .0 00 .9 50 .9 00 .8 5
2

4

6

8

1 0

9 0 n m1 2

1 4

w a v e n u m b e rs [c m

-1

]

1 6

1 8

2 0

2 2

FIG. 2: IR spectra measured during the condensation of evaporated SiO on Si(111). Between two
spectra 12.9 nm of SiO were deposited and the last spectrum corresponds to a final film thicknessof 90 nm. The data was first published in Klevenz et. al. 20105 .

ness increases. The strongest peak can be assigned to a vibration of the oxygen against the
Si atoms, the so-called Si-O strecthing mode (AS1). The second mode at 14.085 µm is connected with the Si-Si strecthing (SS) mode. Those spectra were successfully described using
a dielectric function model according to equation (2) and the derived dielectric function
shows a good accordance to existing literature data based on ellipsometry measurements.5
The development of the relative transmittance with thickness of Mg and Fe is shown in
Fig. 2, respectively. In both cases the expected Drude-like behaviour at higher thickness
can be observed as soon as the percolation is reached after an island-like growth in the
beginning. We also deposited a Fe film at 110 K and were able to reproduce the results byFederov et. al.7 that demonstrated a strong quantum size effect in ultrathin Fe layers that
leads to an increased IR absorbance.
First studies on multilayer structures show a Fano-type interaction of the metal films withthe main vibrational band of SiO that leads to an additional feature in the 10 µm region.
Detailed analysis of this effect and further measurements on more complex structures arework in progress.

4r e l. tr a n s m itta n c e

33

1 .0

1 .0

0 .9

0 .9

0 .8

0 .8

0 .7

0 .7

0 .6

0 .6

0 .5

0 .5

0 .4

0 .4

0 .3

a )

4 .0 n m

M g o n S i( 1 1 1 )

4

6

8

1 0

w a v e l e n g t h [ µm ]1 2

3 .8 n m

0 .3

1 4

b )

4

F e o n S i( 1 1 1 )
6

8

1 0

w a v e l e n g t h [ µm ]

1 2

1 4

FIG. 3: In a) and b) IR spectra measured during the condensation of evaporated Mg and Fe on
Si(111) at room temperature are shown, respectively. The final film thickness of each measurement
is given in the graph. For the evaporation of Mg between two adjacent spectra 0.29 nm were
deposited and in the case of Fe evaporation 0.39nm, respectively.
Acknowledgment

This research has been supported in part by the Deutsche Forschungsgemeinschaft DFG
through grant Schwerpunktprogramm 1385 ”The First 10 Million Years of the Solar System- a Planetary Materials Approach”. The authors acknowledge fruitful discussions with H.-P.
Gail. S. Wetzel and A. Tamanai gratefully acknowledge financial support by T. Yamamotoand by the DFG through grant FOR759 to attend the grain formation workshop.

1

M. Klevenz, S. Wetzel, M. Trieloff, H.-P. Gail, and A. Pucci, phys. stat. sol. (b) 247, 2179 (2010).

2

D. W. Berreman, Phys. Rev. 130, 2193 (1963).

3

A. Lehmann, phys. stat. sol. (b) 148, 401 (1988).

4

U. Teschner and K. Hübner, phys. stat. sol. (b) 159, 917 (1990).

5

M. Klevenz, S. Wetzel, M. Möller, and A. Pucci, Appl. Spectrosc. 64, 298 (2010).

6

R. Brendel and D. Bormann, J. Appl. Phys. 71, 1 (1992).

7

D. V. Fedorov, G. Fahsold, A. Pucci, P. Zahn, and I. Mertig, Phys. Rev. B 75, 245427 (2007).

5

34

Surface structures of presolar alumina grains from unequilibrated
ordinary chondrites
A. Takigawa1, S. Tachibana1, K. Nagashima2, K. Makide2,3, G. R. Huss2, A. N. Krot2, H. Nagahara1,
and K. Ozawa1
1

Department of Earth & Planetary Science, The University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan.

takigawa@eps.s.u-tokyo.ac.jp, 2 Hawai‘i Institute of Geophysics and Planetology, Univerity of Hawai‘i at Mānoa,3

Department of Chemistry, Gakushuin University

1.

Introduction

Corundum (α-Al2O3) is one of the first refractory minerals that condense directly from gas of the solar
composition (Wood & Hashimoto 1993) [1].

Presolar alumina grains, which have highly unusual

oxygen isotopic compositions compared to solar-system materials, are condensates in outflows fromoxygen-rich AGB stars and/or supernovae (e.g., Nittler et al. (1997) [2]).

It is therefore important tounderstand the formation conditions of corundum grains in order to understand the first stage of dust
formation around evolved stars.

The morphological and crystallographic features of presolar alumina

grains should reflect the formation conditions and subsequent thermal history of the grains (Takigawa etal. (2009a, b) [3,4]).

Moreover, morphology of presolar alumina grains may provide clues to understandthe 13-µm features of oxygen-rich evolved stars, which seem to be consistent with disk-like corundum
grains flattened to the crystallographic c-axis (Takigawa et al. (2009a) [3]).

However, morphologicaland crystallographical information of presolar alumina is very limited (Choi et al. 1998; Stroud et al.
2004 [5,6]) because the abundance of presolar alumina in chondrites is low (1−100 ppm) and because
most of grains were found by ion mapping [e.g., Nittler et al. (1997) [2]], which destroys grainmorphology due to ion beam sputtering.

In this study, in order to understand the morphology andcrystallographic orientation of presolar alumina grains, we made detailed observations of morphology andcrystallography of alumina grains from residues unequlibrated ordinary chondrites using field-emission
scanning electron microscopy (FE-SEM) and electron back-scattered diffraction (EBSD) prior to theoxygen isotopic analyses to identify the circumstellar condensates.

2.

Samples and analytical techniques

The acid-residues of ordinary chondrites, Semarkona (LL3.0), Bishunpur (LL3.1), and Roosevelt
County (RC) 075 (LL3.2), prepared by Huss et al. (1995, 2003) [7,8] were used for this study.

The

residues, consisting of acid-resistant grains such as corundum, spinel, hibonite, and SiC, were dispersed

35

onto gold substrates cleaned by ion milling.

Alumina candidate grains were found withcathodoluminescence imaging by a field emission electron microprobe, and 198 corundum grains (0.5−2µm in size) were confirmed by EDS (energy dispersive spectroscopy) equipped to a field emission
scanning microscope, among which 31, 74, and 93 grains were from Semarkona, Bishunpur, and RC 075,
respectively.

For individual grains, secondary electron images were taken from four different directions

and crystallographical information was obtained by EBSD.
Oxygen isotopic compositions of 111 grains were measured with UH Cameca ims-1280 ion
microprobe. The details of analytical technique are described in Makide et al. (2009) [9].

A 0.5−2.5 nA

Cs+ primary ion beam was defocused to ~30 µm for uniform sputtering of grains having different
morphologies and sizes (0.5−3 µm).

Α field aperture of 1500 × 1500 µm2 corresponding to ∼10 µm on

the sample was used to minimize contribution of oxygen signals from the substrate and other grains
surrounding the grain of interest.

Secondary 16O−, 17O−, and 18O− ions were measured in multicollection

mode with the magnetic field controlled by a nuclear magnetic resonance probe.

16

O− and 18O− were

measured by a multicollector Faraday cup and EM (electron multiplier), respectively, with mass resolving
power (MRP) of ∼ 2000, and 17O− was measured using the axial monocollector EM with MRP of ∼ 5600
to separate the interfering 16OH− signal.
>~106 cps.

The primary beam current was adjusted to obtain 16O− signal of

Acquisition time was typically 200 s (10 s × 20 cycles), but it was reduced to 50 – 130 s (10s × 5−13 cycles) during analyses for grains showing presolar signatures of oxygen isotopes to keep grainsfrom being completely sputtered.

After the isotopic measurements, the presolar alumina grains were

reexamined by FE-SEM and EBSD.

3. Results
The 198 alumina grains discovered were classified into three types (A, B, and C) according to theirmorphology.

Type A grains have smooth surfaces or crystal faces (98 grains; Fig. 1a).

Type Bgrains are irregularly shaped and have rough surfaces with 10 to 100 nm-sized fine structures
without crystal faces (67 grains; Fig. 1b).

Type C grains have both type A and type B surface

structures (33 grains), which include two grains with clear bilateral surface structures having crystal
faces on half of the grain surface and irregularly fine structures on another half of the surface (Fig.
1c).

EBSD measurements for 170 grains showed that 87 % of grains are crystalline α-Al2O3(corundum).
Oxygen isotopic compositions of most grains plot along the Carbonaceous Chondrite Anhydrous
Mineral (CCAM) line.

Nine presolar alumina grains were found (~8% of measured grains): sevengrains have oxygen isotopic compositions with positive 17O excesses and small 18O depletions (Group I);

36

two grains show 17O depletions (Group III).
A grains.

The presolar grains consist of 7 type B (Fig. 2) and 2 type

We note that both type A grains were attached to tens-nm-sized type B grains that could be

carriers of presolar isotopic signatures.

The observed presolar/solar corundum ratio of 8% is similar tothat of 7% by Choi et al. (1998) [5], but higher than those of 3% and <0.8 % reported by Nittler et al.
(1997) [2] and Makide et al. (2009) [9], respectively, in spite of essentially the same chemical treatmentof samples and procedure of measurements as ours.
in size of grains analyzed.

This difference could be attributed to the difference

The size of grains measured in our study and in Choi et al. (1998) [5] is about

1 µm on average, but only larger grains (1-5 µm) were measured by Makide et al. (2009) [9].

Most ofpresolar alumina in Tieschitz (H3.6), which were measured by ion-mapping technique, were reported tobe <1.5 µm (Nittler et al. 1997 [2]) being consistent with higher abundances of presolar alumina in ourstudy and in Choi et al. (1998) [5] than that in Makide et al. (2009) [9].

Fig. 1grains.

Alumina grains extracted from unequilibrated ordinary chondrites.

a) type A, b) type B, and c) type C

37

Fig. 2Presolar alumina grain (type B) (RC 075 59-08).

4. Origin of the type B surface structures

The surfaces of type B alumina grains consist of very fine structures (10-100 nm in size).
Such surface features are observed not only for presolar alumina but also for one third of solar
alumina grains.

The type B surface could be the original structure formed in space, but one

may suspect that the acid treatments to extract alumina and other acid-resistant grains from bulk
chondrites alter the surface structures of alumina grains by dissolution.

In order to evaluate thetype B surface structures are original formed in space or artifacts during acid treatments of themeteorites and to assess the resistance of alumina polymorphs to acid treatments, dissolution
experiments of several kinds of crystalline and amorphous alumina particles were conducted.
The dissolution experiments of alumina particles were conducted following the procedure ofacid treatments by Huss & Lewis (1995) [7], based on which the acid residues of chondrites
used in this study were prepared.

In this study, three acid treatments were performed for

α-alumina (corundum); (1) 12M HF - 6M HCl at 25°C for 10 days, (2) 2M H2SO4 + 0.5N
K2Cr7O4 at 75°C for 12 hours, and (3) HClO4 at 190-200°C for 2 hours.

No recognizablechange of surface structures was seen for the α-alumina powders that experienced the acidtreatment (1)-(3) in the order.

This leads to the conclusion that the surface structures of

meteoritic corundum grains in this study are not the artifact formed during acid treatments but
the original features reflecting processes in space.

A possible process that formed the presolar

38

alumina with type B (fluffy and fine-structured) surfaces could be surface sputtering of grains insupernova-induced shockwaves in interstellar medium or the protosolar disk (Jones et al. 1994[10]).

References[1] Wood, J. A. & Hashimoto, A. (1993), “Mineral equilibrium in fractionated nebular systems”,
Geochimica et Cosmochimica Acta, 57, 2377
[2] Nittler, L. R., Alexander, C. M. O'd., Gao, X., Walker, R. M., & Zinner, E. (1997), “Stellar Sapphires:
The Properties and Origins of Presolar Al2O3 in Meteorites”, Astrophysical Journal, 483, 475[3] Takigawa A., Tachibana S., Nagahara H. and Ozawa K. (2009) “Condensation anisotropy ofcorundum around AGB stars and its effect on infrared spectra”, Lunar & Planetary Science, XXXX,
abstract #1731[4] Takigawa, A., Tachibana, S., Nagahara, H. Ozawa, K., and Yokoyama, M. (2009), “Anisotropic
Evaporation of Forsterite and Its Implication for Dust Formation Conditions in Circumstellar
Environments”, Astrophysical Journal, 707, L97[5] Choi, B.-G.-, Huss, G. R., Wasserburg, G. J., and Gallino, R. (1998), “Presolar alumina and Spinel in
Ordinary Chondrites: Origins from AGB Stars and a Supernova”, Science, 282, 1284
[6] Stroud, R. M., Nittler, L. R., & Alexander, C. M. O’D. (2004), “Polymorphism in Presolar Al2O3
Grains from Asymptotic Giant Branch Stars”, Science, 305, 1455
[7] Huss, G. R. & Lewis, R. S. (1995), “Presolar diamond, Sic, and graphite in primitive chondrites:
Abundances as a function of meteorite class and petrologic type”, Geochimica et Cosmochimica Acta,
59, 115[8] Huss, G. R., Meshik, A. P., Smimth, J. B., & Hohenberg, C. M. (2003), “Presolar diamond, siliconcarbide, and graphite in carbonaceous chondrites: Implications for thermal processing in the solar
nebula”, Geochimica et Cosmochimica Acta, 67, 4823
[9] Makide, K., Nagashima, K., Krot, A. N., & Huss, G. R. (2009), “Oxygen Isotopic Compositions ofSolar Corundum Grains”, Astrophysical Journal, 706, 142[11] Jones, A. P., Tielens, A. G. G. M., Hollenbach, D. J., & McKee, C. F. (1994), “Grain destruction inshocks in the interstellar medium”, Astrophysical Journal, 433, 797

39

Grain Formation Workshop 2011

Star formation and dust extinction properties of local galaxies
seen from AKARI
A. Sakurai1 , T.T.Takeuchi1 , F-T.Yuan1 , V. Buat2 , D. Burgarella2
12

Division of Particle and Astrophysical Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan
Laboratoire d’Astrophysique de Marseille, OAMP, Université Aix-Marseille, CNRS, 38 rue Frédéric Joliot-Curie, 13388
Marseille cedex 13, France

Abstract
Accurate estimation of the star formation-related properties of galaxies is crucial forunderstanding the evolution of galaxies. In galaxies, ultraviolet (UV) light emitted
by formed massive stars is attenuated by the dust which is also produced by the SF
activity, and is reemitted at mid- and far- infrared wavelengths (IR). In this study,
we investigated the star formation rate (SFR) and dust extinction using data at UV
and IR. We selected local galaxies which are detected at AKARI FIS 90 µm band. We
measured flux densities at FUV and NUV from the GALEX images. We examined
the SF and extinction of Local galaxies by using four wavebands of AKARI. Then, we
calculated FUV and total IR luminosities, and obtained the SF luminosity (LSF :the totalluminosity related to star formation activity) and the SFR. We found that in most ofgalaxies, LSF is dominated by Ldust . We also found that galaxies with higher SF activity
have a higher fraction of SF hidden by dust. Especially, the SF of galaxies whose SFRs
are > 20 M yr−1 is almost completely hidden by dust. Though these results werequalitatively suggested by previous studies, in this work we could determine the SF andextinction properties with great precision thanks to much larger samples from AKARI
and GALEX all sky surveys.

1 Introduction
Accurate estimation of the star formation -related properties of galaxies is crucial
for understanding the evolution of galaxies. In galaxies, ultraviolet (UV) light emitted
by formed massive stars is attenuated by the dust which is also produced by the SF
activity, and is reemitted at mid- and far- infrared wavelengths (IR). In this study, weinvestigated the star formation rate (SFR) and dust extinction using data at UV andIR.
The total mass of newly formed stars in a galaxy per year is referred to as the starformation rate (SFR). In order to know the star romation (SF) activity, massive stars
are good indicators because they have much shorter lifetime (∼ 106 yr) than the age ofgalaxies and the Universe (order of ∼ 1010 yr). Massive stars (OB stars) are hot andemit ultraviolet light. Then, we can in principle obtain the SFR of galaxies directly
by measuring their UV luminosity. However, as mentioned above, stars produce heavy
elements (or metals ) and release them out through the phenomena during the latest
phase of stellar evolutoin, as planetary nebulae, supernova explosions, and some othermass ejection processes (e.g., Asano et al., 2012) .

1

40These metals are not in the state of gas in the interstellar medium (ISM), but form
tiny solid grains, whose typical size is less than 1 µm. This means that the SF activity
in galaxies is always accompanied by the dust formation, except the very first formationof stars in the Universe. Then, the dust grains are gradually mixed with the interstellar
medium. The UV photons from massive stars are easily absorbed and/or scattered by
dust grains and re-emitted as mid- and far-infrared photons. This is referred to as thedust extinction. Then, we should stress that if we only measure the UV photons frommassive stars in galaxies, the SF activity would be seriously underestimated becausesignificant amount of the energy is obscured by dust (e.g., Kennicutt, 1998) .
So, we investigated star formation -related properties from UV and IR data.

2 Observational data
AKARI is an IR satellite and we used the data of AKARI FIS All-Sky Survey which
has four wavebands at FIR, centered at 65 µm (N60 ), 90 µm (WIDE-S ), 140 µm(WIDE-L), 160 µm (N160 ). IRAS is also an IR satellite and its all sky survey wasbrought with four wavebands before AKARI.
GALEX is a UV satellite which has two wavebands, FUV (Far UV, 1350-1750 Å,
λmean = 1530 Å) and NUV (Near UV, 1750-2750 Å, λmean = 2310 Å) with detection
limits of 19.9 mag and 20.8 mag (Morrissey et al., 2007) . We used the data of GALEX
releases GR4/GR5 all sky survey imaging survey (AIS) in this study.
We started from the AKARI FIS bright source catalog (BSC) v.1 from AKARI allsky survey (Yamamura et al., 2010) to construct our galaxy sample. In order to have asecure sample of galaxies with redshift data, we made a cross match of AKARI sources
with the Imperial IRAS-FSC Redshift Catalogue (IIFSCz), a redshift catalog recently
publishedj (Wang and Rowan-Robinson, 2009) .
After some selection processes, we obtained 6674 galaxies as IR selected sampleandused the GALEX photometry.
We measured FUV and NUV flux densities using GALEX GR4/GR5 images for the6674 galaxies. Then we deleted the wrong images and finaly we obtained available 3891
galaxies.

3 Results
We first obtained the total IR (TIR) luminosity LTIR from AKARI FIS four bands.
LTIR = L2band
AKARI = ∆ν(WIDE-S)Lν (90 µm) + ∆ν(WIDE-L)Lν (140 µm) ,∆ν(WIDE-S) = 1.47 × 1012 [Hz]∆ν(WIDE-L) = 0.831 × 1012 [Hz] .
Figure. 1 shows the relation between the FUV luminosity, LFUV from GALEX photometry and LTIR . Then we caliculate the star formation luminosity LSF using the FUV
luminosity LFUV from GALEX photometry and LTIR . The LSF luminosity is expressedasLSF ≡ LFUV + (1 − η)LTIR2

41

Fig.1 The relation between LFUV and LTIR .The diagonal dotted line represents the caseif LFUV equals LTIR .The blue symbols shows the galaxies which have UV flux densities
below the detection limit of GALEX.

Fig.2 The left panel shows the relation between LFUV and LSF . The right panel shows
the relation between Ldust and LSF .

where η is the fraction of the IR emission produced by dust heated by old stars, which
is not related the current star formation. We adopted 30 % for this fraction (Hirashita
et al., 2003). Figure. 2 shows the contribution of LFUV and Ldust (≡ (1 − η)LTIR ) to
the LSF . With an spectral evolutionary synthesis model Starburst99 (Leitherer et al.,1999), and assuming a constant SFR over 108 yr, the Salpeter IMF ((Salpeter, 1955)
mass range 0.1M –100M ), and the solar metallicity, we obtain the relation between
the SFR and LFUV as (Iglesias-Páramo et al., 2006)log SFRFUV = log LFUV − 9.51

3

42

Fig.3 The left panel shows the contribution of the fraction of SFRFUV to the total SFR.
The right panel shows the relation between LSF and LTIR /LFUV .

The relation between SFR and TIR luminosity islog SFRdust = log LTIR − 9.75 + log(1 − η)We assume that all stellar light is absorbed by dust. Then we can the total SFR
contained FUV and IR lumionsity.
SFR = SFRFUV + SFRdust
The fraction of SFRFUV to the total SFR for the sample galaxies is shown in the left
panel of Fig. 3. There is a sudden drop at SFR > 20 M yr−1 . In the relation between
LSF and the extinction rate, LTIR /LFUV is shown in the right panel in Fig. 3. We
expected high dust extinction rate because of IR selected galaxies.

4 Conclusion
We analyzed star formation-related properties of local galaxies by using AKARI and
GALEX data. Summary and conclusions of this study are as follows:
1. The star formation luminosity, LSF , is dominated by the total infrared luminosity
from dust, LTIR .
2. The contribution of the ultraviolet luminosity, LFUV , has a larger scatter tha thatof the contribution of LTIR .
3. It is difficult to estimate the star formation activity only from the relation between
LSF and LFUV .
4. Galaxies with higher SF activity (SFR > 20 M yr−1 ) have a higher fraction of SF
hidden by dust.
5. This analysis clearly proved that there is a class of glaxies that have small
LFUV /LTIR . These galaxies were dropped on the previous study that used only
bright galaxies in AKARI catalog.

443

From this study, we obtain the similar conclusions with Takeuchi et al. (2010), but withmuch smaller errors. The newly constructed data in this study will provide us with alarge potential of new interesting studies.

References
Asano, R. S., Takeuchi, T. T., Hirashita, H., Inoue, A. K., Dust formation history of galaxies: acritical role of metallicity for the dust mass growth by accreting materials in the interstellar
medium, Earth, Planets, & Space, submitted, 2011.
Hirashita, H., Buat, V., & Inoue, A. K., Star formation rate in galaxies from UV, IR, and Hαestimators, Astronomy and Astrophysics, 410, 83–100, 2003.
Iglesias-Páramo, J., et al., Star Formation in the Nearby Universe: The Ultraviolet and Infrared
Points of View, The Astrophysical Journal Supplement Series, 164, 38–51, 2006.
Kennicutt, R. C., Jr., Star Formation in Galaxies Along the Hubble Sequence. Annual Reviewof Astronomy and Astrophysics, 36, 189–232, 1998.
Leitherer, C., Schaerer, D., Goldader, J. D., et al., Starburst99: Synthesis Models for Galaxies
with Active Star Formation, The Astrophysical Journal Supplement Series, 123, 3–40, 1999.
Morrissey, P., et al., The Calibration and Data Products of GALEX, Astronomy and Astrophysics, 173, 682–697, 2007.
Pollo, A., Rybka, P., & Takeuchi, T. T., Star-galaxy separation by far-infrared color-color
diagrams for the AKARI FIS all-sky survey (bright source catalog version β-1, Astronomy
and Astrophysics, 514, A3, 2010.
Salpeter, E. E., The Luminosity Function and Stellar Evolution., The Astrophysical Journal,
121, 161, 1955.
Takeuchi, T. T., Buat, V., Heinis, S., Giovannoli, E., Yuan, F.-T., Iglesias-Páramo, J., Murata,
K. L., & Burgarella, D., Star formation and dust extinction properties of local galaxies fromthe AKARI-GALEX all-sky surveys . First results from the most secure multiband sample
from the far-ultraviolet to the far-infrared, Astronomy and Astrophysics, 514, A4, 2010.
Wang, L., Rowan-Robinson, M., The Imperial IRAS-FSC Redshift Catalogue, Monthly Noticesof the Royal Astronomical Society, 398, 109-118, 2009.
Yamamura I., Makiuti S., Ikeda N., Fukuda Y., Oyabu S., Koga T., White G. J., AKARI/FIS
All-Sky Survey Bright Source Catalogue Version 1.0 Release Note, ISAS/JAXA, 2010.

5

44

A Study on the Star Formation History of Galaxies Revealed by
the UV-IR Bivariate Luminosity Function
Tsutomu T. Takeuchi1 , Akane Sakurai1 , Fang-Ting, Yuan1
Véronique Buat2 , and Denis Burgarella2
1

Division of Particle and Astrophysical Science, Nagoya University, Furo-cho, Chikusa-ku,
Nagoya 464-8602, Japan
2
Laboratoire d’Astrophysique de Marseille, OAMP, Université Aix-Marseille, CNRS, 38 rueFrédéric Joliot-Curie, 13388 Marseille cedex 13, France
Abstract
Star formation in galaxies is tightly connected to the dust grain formation through the production of heavy
elements. Observationally, it is the most direct to measure the star formation rate by the amount of ionizing UVphotons, but dust formed by the star formation scatter and absorb these photons, and finally re-emit the energy at
the IR. Hence, it is necessary to observe galaxies at both UV and IR to explore their star formation in the cosmic
history. In this work, first we introduce a statistical method to handle the two datasets which have been treatedindependently at each wavelength in a unified manner, and reveal the star formation in galaxies hidden by dust.

1

Introduction

Exploring the star formation history of galaxies is one of the most important topics in modern observational cosmology. Especially, the “true” absolute value of the cosmic star formation rate (hereafter SFR) has been of a central
importance to understand the formation and evolution of galaxies.
However, it has been a difficult task for a long time because of dust extinction. Active star formation (SF) isalways accompanied by dust production through various dust grain formation processes related to the final stageof stellar evolution (e.g., Dwek & Scalo, 1980; Dwek, 1998; Nozawa et al., 2003; Takeuchi et al., 2005c; Asanoet al., 2011). Observationally, SFR of galaxies are, in principle, measured by the ultraviolet (UV) luminosity frommassive stars because of their short lifetime (∼ 106 yr) compared with the age of galaxies or the Universe. However, since the UV photons are easily scattered and absorbed by dust grains, SFR of galaxies is always inevitably
affected by dust produced by their own SF activity. Since the absorbed energy is re-emitted at wavelengths of
far-infrared (FIR), we need observations both at UV and FIR to have an unbiased view of their SF (e.g., Buat etal., 2005; Seibert et al., 2005; Cortese et al., 2006; Takeuchi et al., 2005a, 2010; Haines et al., 2011; Bothwell etal., 2011; Hao et al., 2011).
After great effort of many researchers, the cosmic history of the SFR density is gradually converging at 0 <z < 1. This “latter half” of the cosmic SFR is characterized by the rapid decline of the total SFR, especially thedecrease of the contribution of dusty IR galaxies toward z = 0: While at z = 0, the contribution of the SFR
hidden by dust is 50–60 %, it becomes > 90 % at z = 1 (Takeuchi et al., 2005a). This difference of decrease in
SFR obtained from FUV and FIR has been already recognized mainly in IR studies (e.g., Takeuchi et al., 2000,
2001a,b). Later works confirmed this dusty era of the Universe , and revealed that the dominance of the hidden
SF continues even toward higher redshifts (z ∼ 3) (e.g., Murphy et al., 2011; Cucciati et al., 2011).
Then, a natural question arises: what does the different evolution at different wavelengths represent? To address
this problem, it is very important to understand how we select sample galaxies and what we see in them. Each time
we find some relation between different properties, we must understand clearly which is real (physical) and which
is simply due to a selection effect. Structure of the sample is schematically described by the Venn diagram in
Fig. 1. Some are detected at both bands, but some are detected only at one of the observed wavelength and appear
as upper limits at the other wavelength. In previous studies it was often found that various claims were inconsistent
with each other, mainly because they did not construct a well-controled sample of FUV and FIR selected galaxies.
1

45

Figure 1: A Venn diagram describing the structure of UV- and FIR-wavelength observational dataset. Some aredetected at both bands, but some are detected only at one of the observed wavelength and appear as upper limits atthe other wavelength.
Recently, thanks to new large surveys, some attempts to explore the SFR distribution of galaxies in bivariate way
have been made, through the “SFR function” (e.g., Buat et al., 2007, 2009; Haines et al., 2011; Bothwell et al.,2011). These works are based on the FUV and FIR LFs and their sum, but have not yet address their dependence oneach other. To explore the bivariate properties of the SF in galaxies further, the proper tool is the UV-IR bivariate
luminosity function (BLF).
However, constructing a BLF from two-wavelength data is not a trivial task. When we have a complete fluxlimited1 multiwavelength dataset, we can estimate a univariate luminosity function (LF) at each band, but what wewant to know is the dependence structure between luminosities at different bands. Mathematically, this problem isrephrased as follows: can we (re)construct a multivariate probability density function (PDF) satisfying prescribedmarginals? Although there is an infinite number of degrees of freedom to choose the original PDF, if we can modelthe dependence between variables, we can construct such a bivariate PDF. A statistical tool for this problem is theso-called “copula” (see Sec. 2 for the definition).
A copula has been extensively used in financial engineering, for instance, but until recently there were very
few applications to astrophysical problems (e.g., Koen, 2009; Benabed et al., 2009; Scherrer et al., 2010; Sato etal., 2010, 2011). Takeuchi (2010) introduced the copula to the estimation problem of a BLF. In this work, we apply
the copula-based BLF estimation to the FUV-FIR two-wavelength dataset from z = 0 to z = 1, using data from
IRAS, AKARI2 , Spitzer3 , and GALEX4 .
The layout of this paper is as follows. We define a copula, especially the Gaussian copula, and formulate thecopula-based BLF in Sec. 2. We then describe our FUV-FIR data in Sec. 3. In Sec. 4, we first formulate thelikelihood function for the BLF estimation. Then we show the results, and discuss the possible interpretation of
the evolution of the FUV-FIR BLF. Section 5 is devoted to summary.
Throughout the paper we will assume ΩM0 = 0.3, ΩΛ0 = 0.7 and H0 = 70 km s−1 Mpc−1 . The luminosities
are defined as νLν and expressed in solar units assuming L¯ = 3.83 × 1033 erg s−1 .

1 Or any other observational condition.
2 URL: http://www.ir.isas.ac.jp/ASTRO-F/index-e.html.
3 URL: http://www.spitzer.caltech.edu/.
4 URL: http://www.galex.caltech.edu/.

2

46

2

The Bivariate Luminosity Function Based on the Copula

2.1

Copula: general definition

First, we briefly introduce a copula. A copula C(u1 , u2 ) is defined as follows:G(x1 , x2 )= C [F1 (x1 ), F2 (x2 )](1)

where F1 (x1 ) and F2 (x2 ) are two univariate marginal cumulative distribution functions (DFs) and G(x1 , x2 ) isa bivariate DF. We note that all bivariate DFs have this form and we can safely apply this method to any kind ofbivariate DF estimation problem (Takeuchi, 2010). In various applications, we usually know the marginal DFs (orequivalently, PDFs) from the data. Then, the problem reduces to a statistical estimation of a set of parameters todetermine the shape of a copula C(u1 , u2 ). In the form of the PDF,
g(x1 , x2 ) =∂ 2 C [F1 (x1 ), F2 (x2 )]f1 (x1 )f2 (x2 ) ≡∂x1 ∂x2

c [F1 (x1 ), F2 (x2 )] f1 (x1 )f2 (x2 ),(2)

where f1 (x1 ) and f2 (x2 ) are the PDFs of F1 (x1 ) and F2 (x2 ), respectively. On the second line, c(u1 , u2 ) is referred
to as a differential copula.

2.2

Gaussian copula

Since the choice of copula is literally unlimited, we have to introduce a guidance principle. In many data analyses
in physics, the most familiar measure of dependence might be the linear correlation coefficient ρ. Mathematicallyspeaking, ρ depends not only on the dependence of two variables but also the marginal distributions, which is
not an ideal property as a dependence measure. Even so, a copula having an explicit dependence on ρ would be
convenient. In this work, we use a copula with this property, the Gaussian copula.
One of the natural candidate with ρ may be a copula related to a bivariate Gaussian DF (for other possibilities,
see Takeuchi, 2010). The Gaussian copula has an explicit dependence on a linear correlation coefficient by its
construction. Let( 2)1x
ψ1 (x) = √ exp −
,
(3)
2
2π
∫ x
Ψ1 =
Ψ(x0 )dx0 ,
(4)
−∞
(
)
1
1
ψ2 (x1 , x2 ; ρ) = √exp − xT Σ−1 x ,2(2π)2 det Σ(5)andΨ2 (x1 , x2 ; ρ) =∫ x1 ∫ x2−∞−∞

ψi (x01 , x02 )dx01 dx02 ,(6)

where x ≡ (x1 , x2 )T , Σ is a covariance matrix
(
Σ≡

)
1 ρ
,
ρ 1

and superscript T stands for the transpose of a matrix or vector.
Then, we define a Gaussian copula C G (u1 , u2 ; ρ) as[]−1
C G (u1 , u2 ; ρ) = Ψ2 Ψ−1
1 (u1 ), Ψ1 (u2 ); ρ .

3

(7)

(8)

47

The density of C G , cG , is obtained as[]∂ 2 Ψ2 Ψ1−1 (u1 ), Ψ−1∂ 2 C G (u1 , u2 ; ρ)1 (u2 ); ρc (u1 , u2 ; ρ) ==∂u1 ∂u2∂u1 ∂u2ψ2 (x1 , x2 ; ρ)=
ψ1 (x1 )ψ1 (x2 )
{
]}
)
1
1[
T (
=√
exp − Ψ−1 Σ−1 − I Ψ−1 ,
2
det Σ
G

(9)

[
]T
where Ψ−1 ≡ Ψ−1 (u1 ), Ψ−1 (u2 ) and I stands for the identity matrix.

2.3

Construction of the FUV-TIR BLF

In this work, we define the luminosity at a certain wavelength band by L ≡ νLν (ν is the corresponding frequency).
Then the LF is defined as a number density of galaxies whose luminosity lies between a logarithmic interval[log L, log L + d log L]:
φ(1) (L) ≡dn,d log L(10)

where we denote log x ≡ log10 x and ln x ≡ loge x. For mathematical simplicity, we define the LF as being
normalized, i.e.,∫φ(1) (L)d log L = 1 .
(11)
Hence, this corresponds to a PDF. We also define a cumulative LF as∫ log L
Φ(1) (L) ≡φ(1) (L0 )d log L0 ,(12)

log Lminwhere Lmin is the minimum luminosity of galaxies considered. This corresponds to the DF. If we denote univariate(1)(1)
LFs as φ1 (L1 ) and φ2 (L2 ), then the BLF φ(2) (L1 , L2 ) is described by a differential copula c(u1 , u2 ) as
]
[
(1)
(1)
(13)
φ(2) (L1 , L2 ) = c φ1 (L1 ), φ2 (L2 ) .
For the Gaussian copula, the BLF is obtained as
φ

(2)

{
}
) −1 ]
1
1 [ −1 T ( −1(1)
(1)
(L1 , L2 ; ρ) = √exp − ΨΣ −I Ψ
φ1 (L1 )φ2 (L2 ) ,
2
det Σ

(14)

where
[
(
)
(
)]T
(1)
(1)
Ψ−1 = Ψ−1 Φ1 (L1 ) , Ψ−1 Φ2 (L2 ).

(15)

For the UV, we adopt the Schechter function (Schechter, 1976).
)1−α1
[ (
)]
(
L
L
(1)
φ1 (L) = (ln 10) φ∗1
exp −
.
L∗1
L∗1
(16)
For the IR, we use the analytic form for the LF proposed by Saunders et al. (1990),
{
)1−α2
[ (
)]2 }
(
L
1
L
(1)
φ2 (L) = φ∗2
exp − 2 log 1 +
.
L∗2
2σ
L∗2

(17)

We use the re-normalized version of eqs. (17) and (16) so that they can be regarded as PDFs, as mentioned above.
4

48

2.4

Selection effects: another benefit of a copula BLF

Another benefit of copula is that it is easy to incorporate observational selection effects which always exist in any
kind of astronomical data. In a bi(multi)variate analysis, there are two categories of observational selection effects.
1. Truncation
We do not know if a source would exist below a detection limit.
2. Censoring
We know there is a source, but we have only an upper (sometimes lower) limit for a certain observable.
As we mentioned above, we have to deal with both of these selection effects carefully to construct a BLF fromdata at the same time. It is terribly difficult to incorporate these effects by heuristic methods, particularly for a
nonparametric methods (Takeuchi 2012, in preparation). In contrast, since we have an analytic form for the BLF,
the treatment of upper limits is much more straightforward (Takeuchi, 2010). We see how it is treated in thelikelihood function in Sec. 4.

3

Data

3.1

FUV-FIR data construction

We have constructed a dataset of galaxies selected at FUV and FIR by GALEX and IRAS for z = 0. At higher
redshifts, GALEX and Spitzer data are used for z = 0.7 and EIS and Spitzer for z = 1.0 in the Chandra Deep
Field South (CDFS). We explain the details of each sample in the following.
In the Local Universe, we used the sample compiled by Buat et al. (2007). This sample was constructedbased on IRAS all-sky survey and GALEX All-Sky Imaging Survey (AIS). This dataset consists of UV- and IRselected samples. The UV-selection was made by GALEX FUV (1530 Å) band with FUV < 17 mag (hereafter,
all magnitudes are presented in AB mag). This corresponds to the luminosity lower limit of LFUV > 108 L¯ .
Redshift information was taken from HyperLEDA (Paturel et al., 2003) and NED. The IR-selection is based on
the IRAS PSCz (Saunders et al., 2000). The detection limit of the PSCz is S60 > 0.6 Jy. Redshifts of all PSCz
galaxies were measured. All UV fluxes were remeasured with the package we have developed (Iglesias-Páramo etal., 2006) to avoid the shredding of galaxies. Details of the sample construction are explained in Buat et al. (2007).
The number of sample galaxies is 644.
We also constructed a new, much larger sample of IR-selected galaxies by AKARI FIS All-Sky Survey. We
started from the AKARI FIS bright source catalog (BSC) v.1 from AKARI all sky survey (Yamamura et al., 2010)
to construct this sample. This sample is selected at WIDE-S band (λ = 90 µm) of the AKARI FIS (Kawada etal., 2007). The detection limit is S90 > 0.2 Jy. We first selected AKARI sources in the SDSS footprints. Then, tohave a secure sample of galaxies with redshift data, we made a cross match of AKARI sources with the Imperial
IRAS-FSC Redshift Catalogue (IIFSCz), a redshift catalog recently published (Wang & Rowan-Robinson, 2009),
with a search radius of 36 arcsec.. Since about 90 % of galaxies in the IIFSCz have spectroscopic or photometricredshifts at S60 > 0.36 Jy, the depth of the sample is defined by this matching. We measured the FUV and NUVflux densities with the same procedure as the IRAS-GALEX sample. The detection limits at FUV and NUV of
this sample are 19.9 mag and 20.8 mag (Morrissey et al., 2007). A corresponding UV-selected sample is underconstruction, hence we have only IR-selected sample. The number of galaxies is 3,891. For more information on
this sample and properties of galaxies, see Sakurai et al. (2011).
At higher-z, our samples are selected in the CDFS. GALEX observed this field at FUV and NUV (2300 Å) as
a part of its deep imaging survey. We restricted the field to the subfield observed by Spitzer/MIPS as a part of the
GOODS key program (Elbaz et al., 2007) to have the corresponding IR data. The area of the region is 0.068 deg2 .
Precise description of our high-z samples are presented in Buat et al. (2009) and Burgarella et al. (2006).
At z = 0.7, the NUV-band corresponds to the FUV rest frame of the sample galaxies. We thus constructed
the sample based on NUV-selection. Redshifts were taken from the COMBO-17 survey (Wolf et al., 2004). Data
reduction and redshift association are explained in Burgarella et al. (2006). We truncated the sample at NUV =25.3 mag so that more than 90 % of the GALEX sources are identified in COMBO-17 with R < 24 mag. We
set the MIPS 24 µm upper limit as 0.025 mJy. For the IR-selected sample, we based on the GOODS sample and5

49

matched the GALEX and COMBO-17 sources. The sizes of both the UV- and IR-selected samples are 340 and
470, respectively.
In contrast to z = 0.7, since NUV-band of GALEX corresponds to 1155 Å in the restframe of galaxies at
z = 1.0, we cannot use NUV as the primary selection band as restframe FUV. Instead, we start from U -band
selection. The CDFS/GOODS field is covered by the EIS survey (Arnouts et al., 2001). We then cross-matched
the U -band sources with the COMBO-1 sample to have redshifts. We set the limit at U = 24.3 mag to avoid
source confusion. The IR flux densities were taken from Spitzer MIPS 24 µm data, and the same upper limit is put
for non-detections. Again the IR-selected sample was constructed from the GOODS, as the z = 0.7 sample. The
sizes of the z = 1.0 UV and IR-selected samples are 319 and 1033, respectively.

3.2

Far-UV and total IR luminosities

We are interested in the SF activity of galaxies and its evolution. Hence, luminosities of galaxies representativeof SF activity would be ideal. As for the directly visible SF, obviously the UV emission is appropriate for thispurpose. We define the FUV luminosity of galaxies LFUV as LFUV ≡ νLν @FUV. For z = 0 galaxies, FUV
corresponds to 1530 Å. At higher redshifts, as we explained, LFUV is calculated from NUV flux density at z = 0.7
and U -band flux density at z = 1.0, respectively.
In contrast, at IR, the luminosity related to the SF activity is the one integrated over a whole range of IR
wavelengths (λ = 8–1000 µm), LTIR . For the Local IRAS-GALEX sample, it is quite straightforward to define
LTIR since the IRAS galaxies are selected at 60 µm. We adopted a formula LTIR = 2.5νLν @60 µm. This roughapproximation is, in fact, justified by SEDs of galaxies with ISO 160 µm observations (Takeuchi et al., 2006).
Also, for AKARI-GALEX sample, we have excellent flux density data at FIR. We used the following TIR
estimation formula from AKARI two wide bands (Takeuchi et al., 2005b, 2010),
L2bandAKARI = ∆ν(WIDE-S)Lν (90 µm) + ∆ν(WIDE-L)Lν (140 µm),

(18)

where
∆ν(WIDE-S) = 1.47 × 1012 [Hz] ,(19)∆ν(WIDE-L) = 0.831 × 1012 [Hz].

However, since deep FIR data of higher-z galaxies are not easily available to date, we should rely on the datataken by Spitzer MIPS 24 µm. we should use conversion formulae from MIR luminosity to LTIR :log LTIR [L¯ ] =log LTIR [L¯ ] =1.23 + 0.972 log L15 [L¯ ] ,(20)
22.27 + 0.707 log L12 [L¯ ] + 0.0140 (log L12 [L¯ ]) ,(21)

which are an updated version of the formulae proposed by Takeuchi et al. (2005b) and also used by Buat et al.
(2009). Here L15 and L12 are luminosities νLν @15µm and 12µm, i.e., 24/(1+z) at z = 0.7 and 1.0, respectively.
The estimated LTIR slightly depends on which kind of conversion formula is used, but for our current purpose,it does not affect our conclusion and we do not discuss extensively here. We can find intercomparison of the
MIR-TIR conversion formula in Buat et al. (2009).
This situation will be greatly improved with Herschel5 data. We will leave this as our future work with Herschel
H-ATLAS (Eales et al., 2010) and GAMA (Driver et al., 2011) data (Takeuchi et al. 2012, in preparation).

4

Results and Discussion

4.1

FUV and TIR univariate LFs

In order to estimate the FUV-TIR BLF, first we have to examine our setting for the FUV and TIR univariate LFs.
The IRAS-GALEX sample, the validity of the Local univariate LFs are already proved (see Fig. 3 of Buat et al.,2007). Hence, we can safely use the analytic formulae of FUV and TIR LFs at z = 0 [eqs. (16) and (17)].
We use the Schechter parameters presented by Wyder et al. (2005) for GALEX FUV: (α1 , L∗1 , φ∗1 ) = (1.21, 1.81×
109 h−2 [L¯ ], 1.35 × 10−2 h3 [Mpc−3 ]). For the TIR, we used the parameters (α2 , L∗2 , φ∗2 , σ) = (1.23, 4.34 ×
5 URL: http://herschel.esac.esa.int/.

6

50

Figure 2: The FUV and TIR LFs at z = 0.7 (left and right panel, respectively). Open squares are a nonparametric
LFs estimated from our CDFS multiwavelength data with GALEX NUV- and Spitzer MIPS 24 µm-selections.
Dotted lines represent the FUV LFs at this redshift bin taken from Arnouts et al. (2005). Dashed lines depict the
TIR LFs derived from the evolutionary parameters at z = 0.7 given by Le Floc’h et al. (2005). Because of awell-known large density enhancement at this redshift, we renormalized to remove the overdensity.
108 h−2 [L¯ ], 2.34 × 10−2 h3 [Mpc−3 ], 0.724) (Takeuchi et al., 2003) obtained from the IRAS PSCz galaxies(Saunders et al., 2000), and multiplied L∗1 with 2.5 to convert the 60-µm LF to the TIR LF.
For higher redshifts, ideally we should estimate from the BLF estimation simultaneously. It is, however, quitedifficult for our current samples because of the limited number of galaxies. We instead used the LF parameters at
z = 0.7 and 1.0 obtained by previous studies on univariate LFs and modeled the FUV and TIR univariate LFs atthese redshifts and examined their validity with nonparametric LFs estimated from the data. We use the parameters
compiled by Takeuchi et al. (2005a). Parameters of the evolution of the TIR LF are obtained by approximating theevolution in the form
[
]
L2
(1)
(1)
φ2 (L2 , z) = g(z)φ2,0
(22)
f(z)
(1)

where φ2,0 (z) is the local functional form of the TIR LF. Le Floc’h et al. (2005) assumed a power-law form for theevolution functions asf(z) = (1 + z)Q , g(z) = (1 + z)P(23)

and obtained P = 0.7 and Q = 3.2, with α remaining constant. The Schechter parameters at z = 0.7 and 1.0 are
directly estimated by Arnouts et al. (2005) and we adopt their values (see Table 1 of Takeuchi et al., 2005a).
Then we estimated the FUV and TIR LFs with the stepwise maximum likelihood method and the variant of
C − method from our CDFS multiwavelength data (for the estimation method, see, e.g., Takeuchi et al., 2000;
Johnston, 2011, and references therein). The obtained univariate LFs at z = 0.7 and 1.0 are presented in Figs. 2
and 3. We also show the analytic model LFs in these figures. We note that a well-known large density enhancement
locates in the CDFS at z = 0.7 (e.g., Salimbeni et al., 2009), we have renormalized the LFs to remove the effect ofoverdensity. In the following analysis, we normalize the univariate LFs according to eq. (11) so that we can treat
the univariate LFs as PDFs, hence this does not affect the following analysis at all.
Because of the small sample size, the LF shape is not perfectly agree with the supposed functional forms, but
the nonparametric LFs acceptably similar to the analytic functions both for the FUV and TIR at each redshift. We
stress that the analytic functions are not the fit to the data, but estimated from other studies. This implies that theestimated evolutionary parameters of the LFs work generally well. Thus, we can use the higher-z univariate LFs
as the marginal PDFs for the estimation of the FUV-TIR BLFs.
7

51

Figure 3: The FUV and TIR LFs at z = 1.0. Symbols are the same as in Fig. 2, but FUV samples are selected at
U -band.

4.2

Copula likelihood for the BLF estimation

By using the estimated univariate FUV and TIR LFs as given marginals, we can estimate only one parameter,
the linear correlation ρ by maximizing the likelihood function. The structure of a two-band selected data is
ik
ik
k
k
(LiFUV, jUV, LiTIR, jIR), ik = 1 . . . nk . Here, jband (band: FUV or TIR) stands for the upper limit flag: jband = 0:
detection and jband = −1: upper limit. Another index k is the indicator of the selected band, i.e., k = 1 means asample galaxy is selected at UV and k = −1 means it is selected at IR. The likelihood function is as follows.
¯
( k
)
k
¯ ik = 1, · · · nk , k = 1, −1
ln L LiFUV, LiTIR
{
nk
ik
i
∑
∑
[
( k
)](1+jUV)(1+jIRk )k=ln pdet LiFUV, LiTIR
k=

n

1
−1

: UV sel: IR sel

ik =1

ik
ik }
)
)
UV
IR
[ UL: IR ( ik)] (1−k)(−j[ UL: UV ( ik)] (1+k)(−j
ik
ik
2
2
+ ln p
LFUV , LTIR,
+ ln p
LFUV , LTIR

(24)

( k
)
k
where pdet LiFUV, LiTIRis the probability for ik th galaxy to be detected at both bands and to have luminosities
ik
ik
LFUV and LTIR ,
( k
)
k
( ik
)
φ(2) LiFUV, LiTIR
;ρ
ik
det
∫ ∞
,
(25)
p
LFUV , LTIR ≡ ∫ ∞
0
0
0
(2)
0
φ (LFUV , LTIR ; ρ) dLTIR dLFUVLlim
FUV (zik )Llim
TIR (zik )

( k
)
k
k
pUL: UV LiFUV, LiTIRis the probability for ik th galaxy to be detected at IR band and have a luminosity LiTIR, butik
not detected at UV band and only an upper limit LFUV,jk is available,∫ LiFUV,j
k
( k
)
0
k
pUL: UV LiFUV, LiTIR≡ ∫ ∞∫ ∞0

Llim
TIR (zik )

k

(
) 0
k
φ(2) L0FUV , LiTIR
dLFUV

φ(2) (L0FUV , L0TIR ; ρ) dL0TIR dL0FUV

8

,

(26)

52

Figure 4: Schematic BLF. (1) Diagonal: The energy from SF is emitted equally at UV and IR with any SF activity.
(2) Upward: The more active the SF in a galaxy is, the more luminous at the IR (dusty SF). (3) Downward: The
more active the SF is, the more luminous at the UV ( transparent SF).
( k
)
k
k
and pUL: IR LiFUV, LiTIRis the probability for ik th galaxy to be detected at UV band and have a luminosity LiFUV,
ik
but not detected at IR band and only an upper limit LTIR,jk is available,∫ LiTIR,jk( k)φ(2) LiFUV, L0TIR dL0TIR

k

( k
)
k
pUL: IR LiFUV, LiTIR
≡∫ ∞

0∫

Llim
FUV (zik )

∞

φ

(2)

(L0FUV , L0TIR ; ρ) dL0TIR dL0FUV

.

(27)

0

The denominator in eq. (25) is introduced to take into account the truncation in the data by observational flux
selection limits at both bands (e.g., Sandage et al., 1979; Johnston, 2011). We should also note that it often
happens that the same galaxies are included both in UV- and IR-selected sample. In such a case, we should countthe same galaxies only once to avoid double counting of them. Practically, such galaxies are included in any of thesamples, because they are detected at both bands and are symmetric between UV- and IR-selections.

4.3

The BLF and its evolution

Using the Gaussian copula, now we can estimate the BLF. The visible and hidden SFRs should be directly reflectedto this function. Dust is produced by SF activity, but also destroyed by SN blast waves as a result of the SF. Many
physical processes are related to the evolution of the dust amount. Thus, first of all, we should describe statistically
how it evolved, as stated in Introduction.
First, we summarize how to interpret the UV-IR BLF schematically in Fig. 4. First, we see the case thatthe ridge of the BLF is straight and diagonal [see (1) in Fig. 4]. This means that the energy from SF is emitted
equally at UV and IR with any SF activity. If the relation is diagonal but has an offset horizontally or vertically,
this suggests that a constant fraction of energy is absorbed by dust and reemitted. Second, if the ridge is curved
upward, it means that the more active the SF in a galaxy is, the more luminous at the IR [dusty SF: (2) in Fig. 4].

953

Figure 5: The BLF of galaxies from IRAS and GALEX at z = 0. Contours are the analytic model constructed by
a Gaussian copula and univariate FUV and TIR LFs.
Third, if the ridge is bent downward, the more active the SF is, the more luminous at UV [‘transparent’ SF: (3) inFig. 4].
Now we show the estimated BLFs in Figs. 5–8. In the Local Universe, the BLF is quite well constrained. The
estimated correlation coefficient ρ is very high ∼ 0.95, both from IRAS-GALEX and AKARI-GALEX datasets.
The apparent scatter of the LFUV –LTIR is found to be due to the nonlinear shape of the ridge of the BLF. Thisbent shape of the BLF was implied by preceding studies (Martin et al., 2005), and we could quantify this feature.
The copula BLF naturally reproduced it.
At higher redshifts (z = 0.7–1.0), the linear correlation remains tight (ρ ∼ 0.85–0.9), even though it isdifficult to constrain the low-luminosity end from the data in this analysis (Spitzer-GALEX in the CDFS). It will
be interesting to apply this method to better forthcoming data.
Though the whole shape cannot be perfectly determined by the current data, we find that ρ in the copula LF ishigh and remarkably stable with redshifts (from 0.95 at z = 0 to 0.85 at z = 1.0). This implies the evolution of
the UV-IR bivariate LF is mainly due to the different evolution of the univariate LFs, and may not be controled by
the evolution of the dependence structure.
In this study, we did not try to examine if the Gaussian copula would be a proper choice as a model of the
FUV-TIR BLF. At z = 0, for example, when UV-selected data comparably large as the AKARI sample are ready,we will be able to specify (a family of) copulas appropriate for this analysis. At higher redshifts, our current data
are not deep enough to do it. Herschel, ALMA, and SPICA6 data, together with optical deep surveys, will enable
us to explore and constrain the BLF evolution in the whole history of galaxy evolution.

6 URL: http://www.ir.isas.jaxa.jp/SPICA/SPICA HP/index English.html.

10

54

Figure 6: The same as Fig. 5 but from AKARI and GALEX.

5Conclusion

To understand the visible and hidden SF history in the Universe, it is crucial to analyze multiwavelength data in a
unified manner. The copula is an ideal tool to combine two marginal univariate LFs to construct a bivariate LFs. It
is straightforward to extend this method to multivariate DFs.
1. The Gaussian copula LF is sensitive to the linear correlation parameter ρ.
2. Even so, ρ in the copula LF is remarkably stable with redshifts (from 0.95 at z = 0 to 0.85 at z = 1.0).
3. The nonlinear structure of the BLF is naturally reproduced by the Gaussian copula.
4. This implies the evolution of the UV-IR BLF is mainly due to the different evolution of the univariate LFs,
and may not be controled by the dependence structure.
The data from Herschel, ALMA, and SPICA data will improve the estimates drastically, and we expect to specify
the full evolution of the UV-IR BLF in the Universe. We stress that the copula will be a useful tool for any otherkind of bi- (multi-) variate statistical analysis.

Acknowledgments
This work is based on observations with AKARI, a JAXA project with the participation of ESA. This research has
made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory,
California Institute of Technology, under contract with the National Aeronautics and Space Administration. TTT
has been supported by Program for Improvement of Research Environment for Young Researchers from Special

11

55

Figure 7: The same as Fig. 5 but from Spitzer and GALEX at z = 0.7.
Coordination Funds for Promoting Science and Technology, and the Grant-in-Aid for the Scientific Research Fund
(20740105, 23340046) commissioned by the MEXT. TTT, AS, and FTY have been partially supported from theGrand-in-Aid for the Global COE Program “Quest for Fundamental Principles in the Universe: from Particles
to the Solar System and the Cosmos” from the Ministry of Education, Culture, Sports, Science and Technology(MEXT) of Japan. &References
Arnouts, S., Schiminovich, D., Ilbert, O., et al. 2001, A&A, 379, 740
Arnouts, S., Vandame, B., Benoist, C., et al. 2005, ApJ, 619, L43
Asano, R. S., Takeuchi, T. T., Hirashita, H., Inoue, A. K., 2011, Earth, Planets, & Space, submitted
Benabed, K., Cardoso, J.-F., Prunet, S., Hivon, E., 2009, MNRAS, 400, 219
Bothwell, M. S., Kennicutt, R. C., Johnson, B. D., Wu, Y., Lee, J. C., Dale, D., Engelbracht, C., Calzetti, & D.,
Skillman, E., 2011, MNRAS, 415, 1815
Buat, V., & Burgarella, D. 1998, A&A, 334, 772
Buat, V., Iglesias-Páramo, J., Seibert, M., et al. 2005, ApJ, 619, L51
Buat, V., Takeuchi, T. T., Iglesias-Páramo, J., et al. 2007, ApJS, 173, 404
Buat, V., Takeuchi, T. T., Burgarella, D., Giovannoli, E., & Murata, K. L. 2009, A&A, 507, 693

12

56

Figure 8: The same as Fig. 5 but from Spitzer and GALEX at z = 1.0.
Burgarella, D., et al., 2006, A&A, 450, 69
Cortese, L., Boselli, A., Buat, V., Gavazzi, G., Boissier, S., Gil de Paz, A., Seibert, M., Madore, B. F., & Martin,
D. C. 2006, ApJ, 637, 242
Cucciati, O., et al., 2011, A&A, submitted (arXive: astro-ph/1109.1005)
Dale, D. A., Helou, G., Contursi, A., Silbermann, N. A., & Kolhatkar, S., 2001, ApJ, 549, 215
Driver, S. P., et al., 2011, MNRAS, 413, 971
Dwek, E., & Scalo, J. M., 1980,ApJ, 239, 193
Dwek, E., 1998, ApJ, 501, 643
Eales, S., et al., 2010, PASP, 122, 499
Elbaz, D., et al., 2007, A&A, 468, 33
Haines, C. P., Busarello, G., Merluzzi, P., Smith, R. J., Raychaudhury, S., Mercurio, A., & Smith, G. 2011,
MNRAS, 412, 127
Hao, C.-N., Kennicutt, R. C., Johnson, B. D., Calzetti, D., Dale, D. A., & Moustakas, J., D. 2011, ApJ, 741, 124
Iglesias-Páramo, J., et al., 2006, ApJS, 164, 38
Johnston, R., 2011, A&A Review, 19, 41
Kawada, M., et al., 2007, PASJ, 59, 389
13

57

Koen, C., 2009, MNRAS, 393, 1370
Le Floc’h, E., et al., 2005, ApJ, 632, 169
Martin, D. C., et al. 2005, ApJ, 619, L59
Morrissey, P., et al., 2007, ApJS, 173, 682
Murphy, E. J., Chary, R.-R., Dickinson, M., Pope, A., & Frayer, D. T., 2011, ApJ, 732, 126
Nozawa, T., Kozasa, T., Umeda, H., Maeda, K., & Nomoto, K., 2003,ApJ, 598, 785
Paturel, G., Petit, C., Prugniel, P., Theureau, G., Rousseau, J., Brouty, M., Dubois, P., & Cambrésy, L., 2003,
A&A, 412, 45
Sakurai, A., Takeuchi, T. T., Yuan, F.-T., Buat, V., & Burgarella, D., 2011, Earth, Planets & Space, submitted
Salimbeni, S., et al., 2009, A&A, 501, 865
Sandage, A., Tammann, G. A., & Yahil, A., 1979,ApJ, 232, 352
Sato, M., Ichiki, K., & Takeuchi, T. T., 2010, PRL, 105, 251301
Sato, M., Ichiki, K., & Takeuchi, T. T., 2011, PRD, 83, 023501
Saunders, W., Rowan-Robinson, M., Lawrence, A., Efstathiou, G., Kaiser, N., Ellis, R. S., & Frenk, C. S., 1990,
MNRAS, 242, 318
Saunders, W., et al., 2000, MNRAS, 317, 55
Schechter, P., 1976, ApJ, 203, 297
Scherrer, R. J., Berlind, A. A., Mao, Q., & McBride, C. K., 2010, ApJ, 708, L9
Seibert, M., et al., 2005, ApJ, 619, L55
Takeuchi, T. T., Yoshikawa, K., & Ishii, T. T. 2000, ApJS, 129, 1
Takeuchi, T. T., Ishii, T. T., Hirashita, H., et al. 2001, PASJ, 53, 37
Takeuchi, T. T., Kawabe, R., Kohno, K., et al. 2001, PASP, 113, 586
Takeuchi, T. T., Yoshikawa, K., & Ishii, T. T. 2003, ApJ, 587, L89
Takeuchi, T. T., Buat, V., & Burgarella, D. 2005, A&A, 440, L17
Takeuchi, T. T., Buat, V., Iglesias-Páramo, J., Boselli, A., & Burgarella, D. 2005, A&A, 432, 423
Takeuchi, T. T., Ishii, T. T., Nozawa, T., Kozasa, T., & Hirashita, H. 2005, MNRAS, 362, 592
Takeuchi, T. T., Ishii, T. T., Dole, H., et al. 2006, A&A, 448, 525
Takeuchi, T. T. 2010, MNRAS, 406, 1830
Takeuchi, T. T., Buat, V., Heinis, S., et al. 2010, A&A, 514, A4Wang, L., & Rowan-Robinson, M., 2009, MNRAS, 398, 109
Wolf, C., et al., 2004, A&A, 421, 913
Wyder, T. K., et al., 2005, ApJ, 619, L15
Yamamura I., Makiuti S., Ikeda N., Fukuda Y., Oyabu S., Koga T., & White G. J., 2010, AKARI/FIS All-Sky
Survey Bright Source Catalogue Version 1.0 Release Note, ISAS/JAXA

14

58

An analytic model for the distributionof matter and dust around galaxies
Shogo Masaki
Department of physics, Nagoya University, shogo.masaki@nagoya-u.jp

ABSTRACTWe develop an analytic halo model for the distribution of matter and dust around galaxies. The model results
are compared with the observed surface matter and dust density profiles measured through gravitational
lensing magnification and reddening of background quasars in Sloan Digital Sky Survey reported by Menard
et al.(2010). We show that our model reproduces, for the first time, both the observed matter and dust
distribution profile very well over a wide range of radial distance of 10- 104 h-1kpc. Through the model
fitting of the mass profile, we determine the typical host halo mass of the sample galaxies used in themeasurement to be 2×1012 h-1Msun. We assume that the dust distribution around a galaxy is describedsimilarly to the mass distribution, but with a sharp cut-off at αRvir where Rvir is the galaxy's virial radius andα is a model parameter. Our model reproduces the observed dust profile well if α is greater than unity, andthus suggests that dust is distributed to over 100kpc - 1Mpc from the galaxies.
1. METHOD
We describe our formulation to calculate the mean surface density profile of matter and dust around galaxies.
Our model is based on the so-called halo model approach (Seljak 2000; Cooray & Sheth 2002). The
approach is known as a powerful method to describe statistical quantities, e.g., the correlation function ofgalaxy and the cross-correlation function of galaxy-matter. In this model, the mean surface density profile isdivided into two terms:
2h
ΣX (R) = Σ1h
X (R) + ΣX (R)where “X” denotes matter or dust and R is the distance in the projected two-dimensional plane. The one-haloterm ΣX1h(R) arises from the central halo, and the two-halo term ΣX2h(R) from the neighbouring halos.
The contribution from an individual galaxy halo with mass Mh to the one-halo term ΣX1h(R) is given by the
projection of the radial halo density profile ρX(r|Mh) along the line-of-sight χΣX (R|Mh ) =�dχ ρX (r =�
χ2 + R2 |Mh ).

The one-halo term is a number-weighted average of ΣX(R|Mh)

Σ1h
X (R) =

1
nhalo

� ∞

dn
ΣX (R|Mh ),
dMh
dMh
Mmin

nhalo =

� ∞

Mmin

dMh

dn
,
dMh

where dn/dMh is the halo mass function and Mmin is the threshold halo mass for the sample galaxies. The
threshold mass corresponds to the typical host halo mass of the observed galaxies.
On the other hand, for the two-halo term, we calculate the halo-matter/dust cross-power spectrum. We
obtain the mean surface density profile from the cross-correlation function converted from the cross-power
spectrum. For the formulation, see e.g., Seljak (2000) and Cooray & Sheth (2002).
We adopt the NFW profile derived from cosmological N-body simulations (Navarro et al. 1997) for theradial mass profileρm (r|Mh ) =ρs.
(r/rs )(1 + r/rs )2

We assume that the dust distribution around a galaxy is described similarly to the mass distribution, but witha sharp cut-off at αRvir59

ρd (r|Mh ) ∝1exp(−r/αRvir )[1 − exp(−r/Rdisk )](r/rs )(1 + r/rs )2

where Rvir is the galaxy's virial radius, Rdisk is size of galactic disk and α is a model parameter. The amplitudeis determined via

� ∞0

dr 4πr2 ρd (r|Mh ) = Md = Γ × Mh

where Md is the total halo dust mass associated with a halo with mass Mh and Γ is the dust-halo mass ratio.
For the dust-halo mass ratio Γ, we propose two models. The first one is constant model, i.e., Γ isindependent of halo mass. The dust-halo mass ratio is simply the global density ratio (Komatsu et al. 2011;
Fukugtia 2011)

Γ = 1.73 × 10−5 = Ωgalaxy dust /Ωm .

Because the heavy elements that constitute dust are produced by stars, it may be reasonable to expect that the
dust mass is proportional to the steller mass. In our second model, we model the halo mass dependence of Γ.
We call the model as mass dependent model. We combine three observed relations; the stellar-halo mass
(Leauthaud et al. 2011), the stellar mass-metallicity (Tremonti et al. 2004) and the stellar-hydrogen mass
(Evoli et al. 2011). The derived Γ is shown and compared with that of the constant model in Fig. 1.

2. RESULTS
We fit our matter distribution model to the observed surface matter density through least chi-squared
minimization. We choose Mmin as a free parameter for this fit. The obtained best-fit threshold mass is

Mmin = (2.1 ± 0.4) × 1012 h−1 M⊙ (1σ),
which is consistent with the value derived by Masaki, Fukugita & Yoshida (2012) using cosmological Nbody simulations. Fig. 2 compares the best-fit model profile with the observed one. The data points witherror bars are taken from Menard et al. (2010, hereafter MSFR). The dotted and dashed lines represent the
one-halo and the two-halo terms respectively. The solid line is the resulting matter profile. Our model
reproduces the observed matter profile remarkably well from ∼10kpc to ∼10Mpc. The obtained value of Mminroughly corresponds to the halo mass of L∗-galaxies. For more detailed discussion on Mmin, see Masaki,
Fukugita & Yoshida (2012).
It is encouraging that our simple physical model for the matter distribution yields an excellent fit to theobservation. We follow the same halo model approach to predict the mean surface dust density profile.
For dust distribution, we have two physical parameters, Mmin and α. We have found that poor constraints
are obtained on the parameters when both of them are treated as free parameters. Because Mmin is alreadyestimated using the observed matter distribution, it is sensible to fit our dust distribution model by treating
only α as a free param- eter. After the same least chi-squared minimization procedure, the obtained best-fit α
for the constant and the mass dependent models are, respectively,

α

=

α

=

1.88+0.693−0.443 (1σ) for constant model,

5.67+2.23−1.44 (1σ) for mass dependent model.

Figs.3 and 4 shows the best-fit dust profile of the constant model with α = 1.88 and that of the massdependent model with α = 5.67. The data points are from MSFR. Both models for Γ reproduce the observed
result fairly well.
Overall, our simple models reproduce the observed dust profile very well. Intriguingly, both our models
suggest α ∼ O(1), i.e., halo dust is distributed over ∼ 100kpc from the galaxy. It is also important that theobserved power law surface density Σd ∝ R−0.8 at R = 10kpc − 10Mpc can be explained with α ∼ O(1). The
apparent large-scale dust distribution is explained by the two-halo contributions.

60

REFERENCES
Cooray, A., & Sheth, R. 2002, PhR, 372, 1
Evoli, C., et al. 2011, ApJ, 743, 45
Fukugita, M. 2011, arXiv:1103.4191
Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJs,192, 18
Leauthaud, A., Tinker, J., Bundy, K., et al. 2011, arXiv:1104.0928
Masaki, S., Fukugita, M. & Yoshida, N. 2012, ApJ, 746, 38
Menard, B., Scranton, R., Fukugita, M., & Richards, G. 2010, MNRAS, 405, 1025 (MSFR)
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
Seljak, U. 2000, MNRAS, 318, 203
Tremonti, C., et al. 2004, ApJ, 613, 898

Figure 1

Figure 2

Figure 3

Figure 461

Evolution of grain size distribution in galaxies
Ryosuke S. Asano1 , Tsutomu T. Takeuchi1 , Hiroyuki Hirashita2 , Akio K. Inoue3
and Takaya Nozawa4
1

Depertment of Particle and Astrophysical Science, Nagoya University, Japan
2
Institute of Astronomy and Astrophysics, Academia Sinica, Taiwan
3College of General Education, Osaka Sangyo University, Japan
4
Institute for the Physics and Mathematics of the Universe, University of Tokyo, Japan

Abstract
We investigate the evolution of the grain size distribution in galaxies by constructing a
dust evolution model taking into account grain size distribution. In order to examine theevolution of the grain size distribution, we consider the following two processes: (i) dust
destruction by SN shocks, and (ii) metal accretion onto the surface of grains (referred to as
grain growth). We found that the grains with radii 0.4 − 0.6 µm become dominant on the
dust mass in galaxies after grain growth occurs. Also, the increase of the number of smaller
grains are suppressed by the SN shocks and metal accretion. The grain size distribution ingalaxies depends strongly on the processes in the ISM.

1Introduction

Dust grains inﬂuence various processes which are important to understand the galaxy evolution.
For example, star formation. Star formation is thought to occur in the cold and dense clouds(molecular clouds). Since hydrogen molecules are mainly formed on the surface of dust grains,the molecular formation rate is really higher for the case with dust (e.g., Hirashita & Ferrara,
2002). Furthermore, dust grains and molecules work as a coolant of the interstellar medium(ISM) where stars can be formed. Thus, dust grains drive star formation in galaxies. Also,
dust grains absorb and scatter the stellar light at shorter wavelengths, and they re-emit the farinfrared and submillimeter wavelengths. Therefore dust aﬀects the spectral energy distributionsof galaxies. Consequently, dust is one of the fundamental factor to understand the evolution ofgalaxies (e.g., Yamasawa et al., 2011).
Dust grains are formed by the condensation of heavy elements. Most of these heavy elements(heavier than helium) released by stellar mass loss during stellar evolution or supernovae (SNe)condense into dust grains. However, dust grains are not only formed from stars but also destroyedby SNe blast waves (e.g., Nozawa et al., 2003; Zhukovska, Gail, & Trieloﬀ, 2008). In addition,
we should consider the grain growth in the ISM by accretion of heavy elements (e.g., Draine,2009). Thus, dust evolution passes through various processes closely related to each other in the
ISM.
Up to now, dust evolution has been studied with various dust evolution model (e.g., Inoue,
2003; Zhukovska, Gail, & Trieloﬀ, 2008; Valiante et al., 2009), but most of them are based on
the assumption of a typical size of grains. However, since the eﬀects of dust grains to the galaxy
evolution, which are already mentioned above, depend on the grain size distribution in galaxies,
the grain size distribution is important information to understand the evolution of galaxies.
162

Grain size distribution produced by asymptotic giant branch stars (AGB stars) are thought
to be the log-normal distribution at the peak of ∼ 0.1 µm (e.g., Winters et al., 1997). As for
SNe II, Nozawa et al. (2007) showed that the grain size distributions produced by SNe are biased
to large sizes (≥ 0.01 µm). Small grains are selectively destroyed by the shocks (Nozawa et al.,2006). Also, Hirashita & Kuo (2011) showed that the contribition of the grain growth to thedust evolution depends on the grain size distribution. Thus, since grain size distributions ingalaxies are aﬀected by various processes simultaneously, we investigate the evolution of dust
mass and the grain size distribution in galaxies by constructing a dust evolution model taking
into account these processes.

2

Dust evolution model

In this study, we consider AGB stars, SNe II, dust destruction by SN shocks, and the grain
growth in the ISM as a process which changes the grain size distribution in galaxies. Here we do
not consider the contribution of the SNe Ia. However, for example, Nozawa et al. (2011) provedthat it is ineﬃcient to the total dust budget in galaxies, and there are no observtional support
for SNe Ia. Hence we do not take into account SNe Ia in this study. In this section, we brieﬂy
introduce four processes mentioned above. Throughout this paper, we assume that the shape ofdust grains is spherical for simplicity. So, we havem(a) =

4πa3
ρj ,
3

(1)

where a is the grain size, and ρj is the bulk density of each species j.

2.1

AGB stars

The grain size produced by AGB stars is uncertain to date. Here, we assume log-normal distribution at the peak of ∼ 0.1 µm (Winters et al., 1997) (left panel of Fig. 1). Recently, Yasuda &Kozasa (2012) have examined the size distribution of SiC grains produced by Carbon-rich AGB
stars, and their results also have suggested that the size distribution is roughly log-normal. As
for the AGB dust mass, we adopt the data of Zhukovska, Gail, & Trieloﬀ (2008).

2.2

SNe IIThe data of dust mass and grain size distribution produced by SNe II are taken from Nozawa etal. (2007) (right panel of Fig. 1). Nozawa et al. (2007) examined the dust mass and grain size
distribution processed by sputtering due to the passage of the forward and reverse shocks. They
calculated two cases for the mixing in the helium core, unmixed and mixed model. In our study,
we adopt the results for unmixed model, since the grain size distributions obtained by unmixed
model are consisitent with observation (e.g., Hirashita et al., 2005). The grain size distribution
is like a broken power-law.

2.3

Dust destruction by SN shocks

Dust grains in the ISM become smaller size or are destroyed by sputtering due to the passageof SN shocks. Thus, in order to trace the evolution of grain size distribution, it is neccesary
to consider dust destruction taking into account the grain size. As for the eﬃciency of dust
destruction for each size, we adopt the results of Nozawa et al. (2006).
2

63

Figure 1: Grain size distribution produced by AGB stars (left panel) or SNe II with 20 M
(right panel).
We formulate the dust destruction as follows. The number of grains with radii between a
and a+da, f (a)da, is expressed as00

f (a)da = ξ(a, a )f (a )da,(2)0where ξ(a, a ) is the conversion eﬃciency, which is deﬁned as the number fraction of grains with000radius between a and a +da eroded to grains with radius between a and a+da by sputtering.
0
0
We note that if a ≤ a, ξ(a, a ) = 0. Thus, the change in the number of dust grains with radii
between a and a+da after the passage of a SN shock, dN (a), is given by
∫ ∞
0
0
0
dN (a) =ξ(a, a )f (a )da − f (a)da.
(3)0

The conversion eﬃciency ξ depends on the ISM density surrounding the SN nSN , and large nSN
increases the destruction eﬃciency of dust grains (Nozawa et al., 2006; Yamasawa et al., 2011).
In this study, we assume nSN = 1.0 cm−3 .
The swept mass by a SN shock Mswept depends on nSN and metallicity in the ISM. Here,
we adopt the ﬁtting formula derived by Yamasawa et al. (2011),
Mswept /M = 1535n−0.202[(Z/Z ) + 0.039]−0.289 ,SN(4)where Z is the metallicity in the ISM. We set Z = 0.02 (Anders & Grevesse, 1989) throughout
this paper.

2.4

Grain growth in the ISM

In this subsection, we formulate the evolution of the grain size distribution by the accretion ofheavy elements onto the surface of grains (grain growth). Here we only treat the grain growth
as the growth of refractory materials. In fact, volatile also grow, but they evapolate if the clouds
disapper. As for volatile, we will treat it in our future work.
We deﬁne f (a, t)da as the number of grains with radii between a and a+da at time t. If
we consider only the accretion of heavy elements, the following equation holds,∂f (a, t)∂+[f (a, t)ȧ] = 0,
∂t
∂a
3

(5)

64

where ȧ ≡ da/dt is the growth rate of the grain radius. Also, the growth rate of the grain mass
is expressed asdm(a)= πa2 αρcl
(6)
Z hvi,
dt
where α is the sticking coeﬃcient, ρcl
Z is the mass density of gaseous metals, and hvi is the meal
velocity of metals in gas phase. And,cl
ρcl
Z = ρISM Z(1 − δ),(7)where ρcl
ISM is the ISM density, and δ is the dust abundance in the metal mass. From Eq. (1),(6), (7), we obtainαρclda
ISM Z(1 − δ)hvi
ȧ ≡
=
.
(8)
dt
4ρj
From this equation, we note that the growth rate does not depend on the grain size.

2.5

Formulation of dust mass evolution for a grain size

From previous subsections, we can formulate the time evolution of the mass of dust species j
3with radii between a and a+da, Md,j (a, t) = 4π
3 a ρj fj (a, t)da, which is expressed as
dMd,j (a, t)= −dt

Md,j (a, t)
SFR(t) + Yd,j (a, t)
MISM (t)][∫ ∞Mswept4π 3000−γSN (t) Md,j (a, t) −a ρjξj (a, a )fj (a , t)daMISM (t)30∂− η [Md,j (a, t)ȧ],
∂a

(9)
(10)
(11)

where fj (a, t)da is the number of the dust grains of species j with radii between a and a+da
at time t, MISM is the total mass in the ISM (total mass of gas, metals and dust), η is themass fraction of the clouds where the grain growth occurs. Here, SFR(t), Yd,j (t), and γSN (t)are the star formation rate, the total dust mass ejected by stars per unit time, and the SN rate,
n . Thus, they are expressed
respectively. For SFR, we assume the schmidt law; SFR(t) ∝ MISM
asMISM (t),τSF∫ 100 Mmd,j (m, Z(t − τm ), a)φ(m)SFR(t)dm,
Yd,j (t) =SFR(t) =

(12)
(13)

mcut (t)∫ 40 MγSN (t) =φ(m)SFR(t)dm,(14)

mcut (t)>8 M

where τSF is the star formation timescale, τm is the lifetime of a star of mass m, md,j (m, Z, a)is the mass of grains with radii a of a star of mass m and metallicity Z, and mcut (t) is the massof a star with lifetime τm = t. In this study, we assume the case with index n = 1 for simplicity.
Since we assume a closed-box model, which assumes the total baryon mass (Mtot ) is constant, if
t = 0, MISM (0) = Mtot . Also, φ(m) is the initial mass function (IMF), and we adopt the Larson
IMF (Larson, 1998) in the mass range (0.1 – 100) M .

4

65

3

Results and Discussions

In Fig. 2, we show the contributions of various propcesses (stars, dust destruction, and the grain
growth) to the grain size distribution in galaxies. We set τSF = 5 Gyr, Mtot = 1010 M , andη = 1.00. It is seen that the dominant sources of dust are stars (AGB stars and SNe II) at thestage of galaxy evolution. However, since the contribution of SN shocks becomes eﬀective att > 1 Gyr (the dust destruction timescale is deﬁned as about 0.1τSF ), small grains are eﬀectively
destroyed by SN shocks after 1 Gyr. Also, heavy elements accrete onto smaller grains because
small grains govern the grain surface. Thus, the increase of the number of smaller grains aresuppressed by the SN shocks and the grain growth. We ﬁnd that if the contribution of thegrain growth becomes eﬀective, the grains with radius 0.4–0.6 µm become dominant on the dust
mass in galaxies. In the dense regions where grain growth occurs, Weingartner & Draine (2001)
suggested that the grain size distribution are biased to large grains. Therefore, the larger grainsare expected to be dominant on the total dust mass in galaxies after the grain growth occurs.
The grain size distribution in galaxies depends strongly on the processes in the ISM. Other
mechanisms which have potential of changing the grain size distribution (especially, shatteringand coagulation) are left in future work.

Acknowledgments.
RSA has been supported from the Grant-in-Aid for JSPS Research under Grant No. 23-5514.
RSA and TTT have been also partially supported from the Grand-in-Aid for the Global COE
Program “Quest for Fundamental Principles in the Universe: from Particles to the Solar Systemand the Cosmos” from the Ministry of Education, Culture, Sports, Science and Technology(MEXT) of Japan. TTT and AKI have been supported by Program for Improvement of Research
Environment for Young Researchers from Special Coordination Funds for Promoting Scienceand Technology, and the Grant-in-Aid for the Scientiﬁc Research Fund (TTT: 23340046, AKI:
19740108) commissioned by the MEXT. HH is supported by NSC grant 99-2112-M-001-006MY3.

References
Anders, E., & Grevesse, N. 1989, GCA, 53, 197
Draine, B. T. 2009, in Cosmic Dust - Near and Far, ASP Conference Series, in press
Hirashita, H., & Ferrara, A. 2002, MNRAS, 337, 921
Hirashita, H., Nozawa, T., Kozasa, T., Ishii, T. T., & Takeuchi, T. T. 2005, MNRAS, 357, 1077
Hirashita, H. & Kuo, T-M. 2011, MNRAS, 416, 1340
Inoue, A. K. 2003, PASJ, 55, 901
Larson, R. B. 1998, MNRAS, 301, 569
Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425,
Nozawa, T., Kozasa, T., Umeda, H., Maeda, K., & Nomoto, K. 2003, ApJ, 598, 785
Nozawa, T. Kozasa, T., & Habe, A. 2006, ApJ, 648, 435
5

66

Nozawa, T., Kozasa, T., Habe, A., Dwek, E., Umeda, H., Tominaga, N., Maeda, K., & Nomoto,K. 2007, ApJ, 666, 955
Nozawa, T., Maeda, K., Kozasa, T., Tanaka, M., Nomoto, K., & Umeda, H. 2011, ApJ, 736, 45
Valiante, R., Schneider, R., Bianchi, S., Andersen, A. C. 2009, MNRAS, 397, 1661
Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296
Winters, J. M., Fleischer, A. J., Le, Bertre, T., & Sedlmayr, E. 1997, A&A, 326, 305
Yamasawa, D., Habe, A., Kozasa, T., Nozawa, T., Hirashita, H., Umeda, H.,& Nomoto, K. 2011,
ApJ, 735, 44
Yasuda, Y., & Kozasa, T. 2012, ApJ, 745, 159
Zhukovska, S., Gail, H. P., Trieloﬀ, M. 2008, A&A, 479, 453

6

67

Figure 2: Contributions of various processes to the grain size distribution with Mtot = 1010 M ,τSF = 5 Gyr and η = 1.00 (bottom panel only). Upper left: stars (AGB stars and SNe II)
only, upper right: stars and SN shocks (dust destruction), bottom: stars, SN shocks, and thegrain growth. Solid, dotted, and dashed lines represent the galactic age t = 0.1, 1.0, 10 Gyr,respectively. Dot-dashed line represent the MRN distribution (Mathis, Rumpl, & Nordsieck,1977) for reference.

7

68

A new diagnostic method to evaluate growth conditionof rock-forming minerals
Hitoshi Miura1,2 , and Katsuo Tsukamoto1
1

Department of Earth Sciences, Tohoku University, 6-3 Aoba, Aramaki, Aoba-ku, Sendai 980-8578,
Japan
2
Corresponding Author E-mail address: miurah@m.tohoku.ac.jp
Abstract
Chondritic meteorites contain rock-forming minerals in which the chemical compositions
are not uniform. The inhomogeneity is termed as zoning proﬁle, which is considered to reﬂect
the growth condition (cooling rate, closed- or open-system, supercooling, and so forth) when
they were formed in early solar system. We investigated the zoning proﬁle, growth velocity,and mineral morphology formed under a constant cooling rate based on a numerical scheme.
From the numerical simulations for a wide range of the cooling rate, we found some simple
power-law relationships in the zoning proﬁle at very beginning of solidiﬁcation. This results
may give us a very useful diagnostic guideline to speculate the cooling rate from zoning
proﬁles recorded in minerals.

1

Introduction

Rock-forming minerals such as olivine, pyroxene, and anorthite show inhomogeneities in thechemical compositions in certain phenocrysts. This is termed as zoning proﬁles, which areconsidered to have been formed during solidiﬁcation of the phenocrysts and/or in alteration
processes in their parent bodies. Therefore, we may obtain signiﬁcant information about theorigins of these minerals by analyzing how their zoning proﬁles were formed. In this paper,
we investigate the formation process of zoning proﬁle, especially at the very beginning ofsolidiﬁcation.
Let us consider that the solid phase comes in contact with the liquid phase at a planar interface in equilibrium (see Fig. 1). Each phase has its own uniform concentration
determined by the phase diagram. The growth of solid is triggered when the equilibrium isperturbed by changes of temperature, pressure, and chemical composition. Onset of growthcauses rejection of incompatible elements from solid to liquid. If the growth velocity is toolarge to homogenize the incompatible element distribution in liquid phase, these elements areconcentrated in front of the advancing solid-liquid interface (boundary layer). The complex
time-dependent behavior of solute redistribution causes signiﬁcant zoning in solid (initialtransient). The solute distribution in the initial transient was formulated by Tiller et al. [1],then extended by Smith et al. [2]. The Smith’ solution suggested that the uptake of incompatible elements into crystal is enhanced as the growth velocity increases. The enhanced
uptake of incompatible elements plays important roles in trace element partitioning between
phenocrysts and host lava [3], incompatible element concentration in melt inclusion [4], andisotopic fractionation during solidiﬁcation [5].
Although the Smith model assumed that the growth velocity is constant for the entire
solidiﬁcation process, it is not necessarily correct in general [6, 7, 8]. The constant growth velocity is certainly an approximation of the complex growth history of crystals in experiments
and in nature [4]. The eﬀect of non-constant growth velocity on the formation of zoning
proﬁle should be investigated because it will give us greater insight to understand the formation of zoned minerals. In addition, it is also important to consider how to speculate theformation condition of mineral from its zoning proﬁle. Previous researches did not provide
any relationship between the growth velocity and the growth conditions, e.g., cooling rate,bulk composition, supercooling, and so forth [3, 4, 5]. The relationship between the growth1

69

Figure 1: Solute redistribution process at the very beginning of solidiﬁcation. (a) A crystal
comes in contact with liquid at the equilibrium condition. Each concentration of impurity incrystal and liquid is determined by the phase diagram at the given temperature. (b) Onset ofgrowth. Rejection of impurity from crystal to liquid makes boundary layer at the advancing
crystal-liquid interface. (c) Further growth. (d) After complete solidiﬁcation. The crystal showsa signiﬁcant zoning termed as initial transient.

2

70

velocity and these growth conditions is essential to consider the formation environment ofminerals.
The aim of this paper is to investigate the formation of zoning proﬁle during initialtransient without an assumption of constant growth velocity. We consider a case that thetemperature decreases at a constant rate from equilibrium temperature at which crystal
coexists with liquid. We numerically solve the crystal growth, elemental partitioning at
the solid-liquid interface, and elemental diﬀusion in liquid phase consistently to obtain thegrowth history and the zoning proﬁle in solid. We carry out the numerical simulations forwide ranges of growth conditions; cooling rate and bulk composition to show some simple
power-law trends with respect to the cooling rate, which may be useful for diagnostics ofzoning proﬁle in minerals.

2

Methods

2.1

Phase-field method

The classical approach to the modeling of ﬁrst order phase transformations involves tracking
the free boundary that separates the growing phase (e.g., crystal) from the parent phase(e.g., supercooled liquid). A phase-ﬁeld method provides an alternative approach, according
to which a new variable, the phase ﬁeld , is introduced to keep track of the phase, taking
on constant values indicative of each of the bulk phases and making a transition between
these values over a thin transition layer that plays the role of the classically sharp interface.
One of the advantages of the phase-ﬁield method is not to require explicit tracking of thefree boundary. This makes the numerical model simpler than the classical approach. The
detailed review of the phase-ﬁeld method would be found in [9], for example.

2.2

Basic equations

The thermodynamically consistent phase-ﬁeld model for a binary alloy based on the entropy
functional is adopted here [10, 11]. This model has been applied for unidirectional solidiﬁcation of binary alloys [12, 13, 8]. We give a uniform temperature ﬁeld T in the present
simulations, so two equations for phase ﬁeld  and concentration ﬁeld C must be solved.
The governing equations are as follows [11]:
]
[
∂
(1)
= Mφ (1 − 16aχg()) {HA (1 − C) + HB C} + 2φ ∇2  ,
∂t
{ [
]}
∂C
vm
= ∇ · D ∇C −C(1 − C)(HB − HA )∇ ,(2)∂tR
where D, vm , and R are the solute diﬀusivity, the molar volume, and the gas constant. a isthe amplitude of stochastic noise that is introduced only at the crystal-liquid interface, andχ is a random number distributed uniformly from -1 to 1 [10, 8]. We use a = 0.3 in thispaper. The solute diﬀusivity D depends on , so we postulate a form
D = DS + p()(DL − DS ),(3)

where DS and DL are diﬀusivities in solid and liquid phases, respectively. Here, HA and HB
are deﬁned by
(
)
1
1
0
HA,B = p ()LA,B−− WA,B g 0 (),
(4)
TA,B
T
where LA,B and TA,B are the latent heat of crystallization per unit volume and the melting
point of pure materials A and B, respectively. We may allow Mφ to depend on composition,and so postulate a form
Mφ = (1 − C)MA + CMB .
(5)

3

71

The variables of WB , MA , MB , and φ are related to the measurable quantities as follows:√√2√ A,B hA,B
3 2 A,B
2 µA,B TA,B
, MA,B =, 2φ = 6 2,(6)WA,B =
2 hA,B TA,B
12 hA,B LA,B
TA,B
where A,B , µA,B , and hA,B are the interfacial energy between solid and liquid, the kinetic
coeﬃcient, and the interface thickness of pure materials A and B. It should be noted that
there is a relationship between hA and hB as
hB
A TB
=
hA
B TA

(7)

because φ does not depend on concentration (see Eq. 6). Finally, g() is the double-wellpotential and p() is a monotonically increasing function which satisﬁes p(0) = 0 and p(1) =1. They are respectively deﬁned byp() = 3 (10 − 15 + 62 ), g() = 2 (1 − )2 .

2.3

(8)

Input parameters

In this study, we adopt olivine as a model mineral. The olivine is one of the major rockforming minerals with the chemical formula (Mg,Fe)2 SiO4 . The ratio of magnesium and ironvaries between the two endmembers of the solid solution series: forsterite (Mg-endmember)and fayalite (Fe-endmember). The concentration C in Eq. (2) means molar fraction of Fewith respect to the total of Mg and Fe. We can ﬁnd the equilibrium composition at agiven temperature from the phase diagram [14]. Input parameters that we adopt in thispaper are listed in Table 1. The diﬀusion coeﬃcient in solid is set zero to investigate theformation of zoning proﬁle without secondary alteration. The numerical result did not change
signiﬁcantly even if we use non-zero DS (≈ 10 9 cm2 s 1 ) because of the very small diﬀusivity
in comparison with in liquid.

2.4

Numerical procedure

Figure 2(a) shows our initial setting for calculation. Initially, a planar seed crystal is inequilibrium with liquid at a temperature of T0 . We consider one-dimensional geometry
assuming that the crystal-liquid interface keeps its ﬂat shape. We adopt a computational
domain of x = 0 − 200 µm and discretize it by 2000 meshes, so the mesh size is ∆x = 0.1 µm.
The initial position of the interface is set at x = 3 µm. The initial concentrations in solid andliquid are given by equilibrium concentrations at a temperature T0 , which are determinedfrom the phase diagram (see Fig. 2b). Decreasing temperature at a constant rate Ṫ causes
disequilibrium at the interface and triggers solidiﬁcation.
We applied ﬁnite-diﬀerence methods to solve Eqs. (1) and (2). The diﬀusion terms inthese equations are solved by implicit method for Eq. (2) and an explicit Euler scheme forEq. (1). No ﬂux boundary condition is applied at boundaries of the computational domain,i.e., the normal derivatives of  and C are zero at these boundaries.

3

Results

3.1

Initial transientFig. 3 shows calculated concentration distribution in solid and liquid phases during solidiﬁcation for a case of cooling rate Ṫ = 5 × 103 K s 1 and initial temperature T0 = 1900 K.
Initially, the solid comes in contact with the liquid at x = 3 µm with their equilibrium
concentrations, CS0 = 0.170 and CL0 = 0.471, respectively (panel a). As the temperature
decreases at a constant rate, the solid-liquid interface advances rightward and the solid region

4

72

Table 1: Input parameters for olivine. Subscripts: A = forsterite, B = fayalite.
notation
vm
DL
DS
LA,B
A
B
TA
TB
µA,B
hA
hB

meaning
molar volume
diﬀusivity in liquid
diﬀusivity in solid
latent heat of crystallizationinterfacial energy
interfacial energy
melting point
melting point
kinetic coeﬃcient
width of interface
width of interfacevalue
45 cm3 mol−1
10−5 cm2 s−1
0 cm2 s−1
1.3 × 1010 erg cm−3
613 erg cm−2
300 erg cm−2
2163 K
1490 K
0.1 cm s−1 K−1
10 × ∆x(A /B )(TB /TA )hA

Figure 2: (a) Initial setting for calculation. (b) Phase diagram of olivine (data from [14]).

5

73

spreads to the liquid side. Since Fe is rejected from the solid to liquid, a boundary layer withhigh concentration is formed immediately away from the advancing interface (panel b). The
concentration in solid increases with the interface advancement. The concentration in solid
exceeds the initial concentration in liquid CL0 at x ≈ 10 µm (panel c). Finally, the systemsolidiﬁes completely and a zoning proﬁle is formed (panel d).
The zoning proﬁle as shown in Fig. 3(d) is divided into two parts: initial transient andsteady growth region. In the initial transient, the proﬁle changes from CS0 to CL0 with acomplex pattern. In a steady growth region, the proﬁle is almost constant at CL0 , namely,
the solid has the same concentration with that of the initial liquid. This is referred as
an eﬀective partitioning coeﬃcient keﬀ being unity, which is deﬁned by keﬀ = CS /CL0 (itshould be noted that the denominator is not at the interface). However, our one-dimensional
calculation will be perturbed before the steady growth is realized because a planar interface
becomes unstable due to morphological instability [15, 12, 8]. Therefore, we do not consider
the steady growth region in the following.
The initial transient is characterized by its complex concentration proﬁle. We notice thesteep increase of concentration just after the beginning of the solidiﬁcation. We deﬁne awidth of initial transient w within which the concentration in solid increases from CS0 to
Cw , where Cw is given by
1 − k0Cw = C0 −
C0 ,
(9)
e
where e is the base of natural logarithm. Using the initial setting, we obtain Cw = 0.360
with k0 = 0.361 (horizontal dashed line in Fig. 3d). It is found that the concentration insolid increases to Cw at x = 7.06 µm. Subtraction of initial position of the interface gives
the width of w = 4.06 µm.

3.2Growth velocity

Fig. 4 shows the position of crystal-liquid interface xi as a function of time. The input
parameters are the same with Fig. 3. Open circles are the results within w (see §3.1 forthe deﬁnition), and ﬁlled ones are after. The slope between two neighboring circles gives
the growth velocity of crystal. It is found that the growth velocity within w is not constantbecause the slope seems to increase gradually: the average growth velocity is calculated as≈ 70 µm s 1 for the ﬁrst 0.01 s and ≈ 190 µm s 1 for the last 0.01 s within w. Aftert ≈ 0.05 s, the growth velocity increases suddenly, which indicates the beginning of thesteady growth region.
Let us assume that the growth velocity within w increases exponentially with respect totime t as given by the following form,
R(t) = R0 et/t0 ,(10)where R0 is the initial growth velocity and t0 is the time constant. Integration of R(t) overt gives the following expression of the position of interface,xi (t) = R0 t0 (et/t0 − 1) + xi0 ,(11)where xi0 = 3 µm is the initial position. The best ﬁt of calculated results within w with
Eq. (11) gives R0 = 64.1 µm and t0 = 0.0256 s. These parameters gives R = 243 µm s 1 at
the end of the region w (t = 0.34 s). The best ﬁt curve is shown by a solid curve in Fig. 4
and agrees with the phase-ﬁeld calculation within w.

3.3

Dependence on cooling rate

As a result of the phase-ﬁeld calculation, we obtain the following three parameters for eachcalculation; width of initial transient w, the initial growth velocity R0 , and the time constant t0 . Fig. 5 shows the dependences of these parameters on the cooling rate Ṫ . Results

6

74

Figure 3: Calculated concentration distribution for a case of Ṫ = 5×103 K s−1 and T0 = 1900 K.
Crystallized region is colored by gray. (a) Initial setting. The initial position of interface is
xi0 = 3 µm. Initial concentrations are CS0 = 0.170 in solid phase and CL0 = 0.471 in liquid
phase, respectively. (b) After t = 0.0328 s. The interface is advancing rightward. Fe-richboundary layer was developed immediately in front of the advancing interface. Zoning proﬁlewas shown to be formed in solid phase. (c) After t = 0.0476 s. The concentration in solidincreases to the initial concentration in liquid CL0 . (d) After complete solidiﬁcation.

7

75

Figure 4: Position of the crystal-liquid interface xi as a function of time. Open circles are theresult of the phase-ﬁeld calculation during initial transient. The region of initial transient iscolored by gray. The solid curve is the best ﬁt result for the initial transient region. The bestﬁt parameters are R0 = 64.1 µm and t0 = 0.0256 s. Filled circles are calculation results after
initial transient. The dashed curve shows a simple extrapolation of the best ﬁt curve.

8

76with diﬀerent initial temperatures T0 , namely, diﬀerent initial compositions, are plotted by
diﬀerent symbols (see legend). It is found that all of these parameters obey power-law trends
with respect to the cooling rate; w ∝ Ṫ 1/2 , R0 ∝ Ṫ 1/2 , and t0 ∝ Ṫ 1 . The power-law trend
is observed for a wide range of T0 except of T0 = 1500 K in which R0 and t0 seem to deviate
from each power-law trend slightly.
For the same Ṫ , w decreases with the decrease of T0 . When Ṫ = 103 K s 1 , we obtainw = 0.43 µm for T0 = 1500 K and w = 12.6 µm for T0 = 2100 K. The diﬀerence of factorof 30 in w suggests that the diﬀerence of 3 orders of magnitude in Ṫ is required to obtainthe same w. For example, the initial transient of 10 µm in width will be reproduced at acooling rate of Ṫ ≈ 103 K s 1 for a Mg-rich composition (T0 = 2100 K), on the other hand,Ṫ ≈ 1 K s 1 for a Fe-rich composition (T0 = 1500 K).
In contrast of w, the growth velocity does not show the monotonic trend on T0 . For thesame Ṫ , R0 decreases as T0 increases from T0 = 1500 K to T0 = 1700 K, however, backs
to increase as T0 increases from T0 = 1700 K to T0 = 2100 K. Similarly, t0 increases as
T0 increases from T0 = 1500 K to T0 = 1800 K, then backs to decrease as T0 increases
from T0 = 1800 K to T0 = 2100 K. The higher growth velocity at large- and small-T0
may be interpreted as follows: the large- and small-T0 mean less impurity, namely, pure
forsterite (Mg2 SiO4 ) or pure fayalite (Fe2 SiO4 ), respectively, so the advancing of interface is
not suppressed by the solute diﬀusion in liquid phase.

4

Conclusion

We numerically simulated solidiﬁcation of olivine system to investigate the formation processof zoning proﬁle at the very beginning of solidiﬁcation, where the iron content in olivine
crystal shows a strong zoning within a certain distance. The strong zoning is termed as
initial transient. When the solidiﬁcation occurs under a constant cooling rate Ṫ , we foundthat the distance of initial transient is proportional to Ṫ 1/2 . During the initial transient,
the growth velocity is not constant but increases exponentially with respect to the time. The
simulation results of growth velocity are ﬁtted with a form of R(t) = R0 et/t0 , where R0 is
the initial growth velocity and t0 is the time constant. We also found that R0 and t0 show
power-law trends with the cooling rate, respectively. The dependences of these parameters
on bulk iron content were also obtained. The power-law trends shown in the initial transient
may give us a very useful diagnostic guideline to speculate the cooling rate from zoning
proﬁles of natural olivines in Earth rocks or chondritic meteorites.

Acknowledgement
We are grateful to Prof. T. Yamamoto for useful discussion on analytic formulations of ournumerical results. H.M. was supported by Tohoku University Global COE Program “Global
Education and Research Center for Earth and Planetary Dynamics”, and by Japan Society
for the Promotion of Science (JSPS), Grant-in-Aid for Young Scientists (B) (23740330).

References[1] W.A Tiller, K.A Jackson, J.W Rutter, and B Chalmers. The redistribution of soluteatoms during the solidiﬁcation of metals. Acta Metallurgica, 1(4):428 – 437, 1953.
[2] V. G. Smith, W. A. Tiller, and J. W. Rutter. A mathematical analysis of solute redistribution during solidiﬁcation. Can. J. Phys., 33:723–745, 1955.
[3] F. Albarede and Y. Bottinga. Kinetic disequilibrium in trace element partitioning between phenocrysts and host lava. Geochimica et Cosmochimica Acta, 36(2):141 – 156,
1972.

9

77

Figure 5: Dependences of initial transient parameters on cooling rate Ṫ ; (a) width of initialtransient w, (b) initial growth velocity R0 , and (c) time constant t0 . The legend indicates thevalue of initial temperature T0 , which stands for the diﬀerence of initial composition of liquid.

10

78

[4] Don Baker. The ﬁdelity of melt inclusions as records of melt composition. Contributions
to Mineralogy and Petrology, 156:377–395, 2008. 10.1007/s00410-008-0291-3.
[5] E. Bruce Watson and Thomas Mller. Non-equilibrium isotopic and elemental fractionation during diﬀusion-controlled crystal growth under static and dynamic conditions.
Chemical Geology, 267(3-4):111 – 124, 2009.
[6] Timothy P. Loomis. An investigation of disequilibrium growth processes of plagioclasein the system anorthite-albite-water by methods of numerical simulation. Contributions
to Mineralogy and Petrology, 76:196–205, 1981. 10.1007/BF00371959.
[7] I. L’Heureux. Oscillatory zoning in crystal growth: A constitutional undercooling mechanism. Physical Review E, 48(6):4460–4469, 1993.
[8] Tomohiro Takaki, Toshimichi Fukuoka, and Yoshihiro Tomita. Phase-ﬁeld simulation
during directional solidiﬁcation of a binary alloy using adaptive ﬁnite element method.
J. Crys. Growth, 283:263–278, 2005.
[9] Robert F. Sekerka. Phase field modeling of crystal growth morphology, pages 176–190.
Perspectives on Inorganic, Organic, and Biological Crystal Growth: From Fundamentals
to Applications. American Institute of Physics, Melville, New York, 2007.
[10] J.A. Warren and W.J. Boettinger. Prediction of dendritic growth and microsegregation
patterns in a binary alloy using the phase-ﬁeld method. Acta Metallurgica et Materialia,
43(2):689 – 703, 1995.
[11] Zhiqiang Bi and Robert F. Sekerka. Phase-ﬁeld model of solidiﬁcation of a binary alloy.
Physica A, 261:95–106, 1998.
[12] William J. Boettinger and James A. Warren. Simulation of the cell to plane front
transition during directional solidiﬁcation at high velocity. Journal of Crystal Growth,200(3-4):583 – 591, 1999.
[13] Zhiqiang Bi and Robert F. Sekerka. Phase ﬁeld modeling of shallow cells during directional solidiﬁcation of a binary alloy. J. Crys. Growth, 237-239:138–143, 2002.
[14] N. L. Bowen and J. F. Schairer. American Journal of Science, 29:151–217, 1935.
[15] W. W. Mullins and R. F. Sekerka. Stability of a planar interface during solidiﬁcationof a dilute binary alloy. J. Appl. Phys., 35:444–451, 1964.

11

79

Alignment of Grains by Radiative Torques:
Discovery of Correlation between Interstellar
Polarization and Dust Temperature
Masafumi Matsumura & Youko Kameura ∗Faculty of Education, Kagawa University, Saiwai-Cho 1-1, Takamatsu, Kagawa 760-8522(Received ; accepted )Abstract
We investigate the eﬃciency of interstellar polarization pλ /Aλ , where pλ is thefractional linear polarization and Aλ is extinction, with using the data obtained by
us with the low-dispersion spectropolarimeter HBS, as well as those in literature.
It is found that the polarization eﬃciency pλ /Aλ is proportional to exp(−β/λ) inwavelength λ ≈ 0.4 − 0.8µm, where β is a parameter which varies from 0.5 to 1.2 µm.
We also ﬁnd that β is negatively correlated with the dust temperature deduced from
infrared data by Schlegel et al. (1998), suggesting that the polarization eﬃciency
is higher in short wavelength for higher temperature. According to the alignment
theory by radiative torques (RATs), if the radiation is stronger and bluer, RATs
will make small grains align better, and the polarization eﬃciency will increase in
short wavelength. Our ﬁnding of the correlation between β and the temperature isconsistent with what is expected with the alignment mechanism by RATs. This paper
is essentially the same as Matsumura et al. (2011), but it includes a model of grains.
Key words: alignment mechanism — ISM: dust, extinction — polarization

1.

Introduction

Interstellar polarization, i.e., dichroic extinction observed in distant stars, shows thatgrains are optically anisotropic and aligned, although the mechanism of the alignment has
been still on debate over more than half century (Lazarian 2007, for a review). The alignment
had been explained with the paramagnetic relaxation of thermally spinning grains that obtain
angular momentum by collisions with gas particles (Davis & Greenstein 1951, hereafter DG).
However, the DG mechanism is not eﬃcient, and it cannot explain the interstellar polarization∗

Present Address: Akita Junior High School, Magodai-machi 72, Kumamoto 861-5254 and Rikigo JuniorHigh School, Sima-machi 5-8-1, Kumamoto 861-4133

1

80

quantitatively.
For more eﬃcient alignment, Purcell (1979) assumed a spin-up of grain by the ejection ofmolecular hydrogen from grain surface (the ”pinwheel mechanism”). If the site of the ejection is
ﬁxed and exist long enough, the grains could be accelerated and rotate suprathermally, i.e. much
faster than the thermal rotation as was assumed in DG, and they could align well. Lazarian &Draine (1999a), however, noticed the phenomenon of relaxation caused by the reorientation ofspin angular momentum of protons by exchanging angular momentum with grain rotation (the”nuclear relaxation”). Once the grain rotates thermally, it remains thermal by this relaxation,even if the pinwheel mechanism is working, and this is called as ”thermal trapping” (Lazarian& Draine 1999a; Lazarian & Draine 1999b). If most of the grains are thermally trapped, thedegree of alignment is low as expected as the DG mechanism.
Jones & Spitzer (1967) showed that the increase of the imaginary part of magnetic susceptibility of grain material will result in better alignment, when paramagnetic grains containferromagnetic clusters (superparamagnetism). However, Roberge & Lazarian (1999) investigated the DG alignment including the Barnett relaxation and other eﬀects, and they showed
that, even if grains have superparamagnetic inclusions, eﬃcient alignment is achieved only when
the grain temperature is much lower than the gas temperature. Thus this mechanism may work
only in limited astronomical environments.
Dolginov & Mitrofanov (1976) ﬁrst pointed out that irregularly shaped grains thathave ”helicity” can spin up by radiative torques (hereafter RATs). More recently, Draine& Weingartner (1996) estimated the RATs with the discrete dipole approximation method,and showed that RATs are very eﬀective to align grains. Since magnetic moments within rotating grains are induced by the Barnett eﬀect, the grains precess around the magnetic ﬁeld,
and then the direction of alignment is usually parallel to the interstellar magnetic ﬁeld1 (e.g.
Draine & Weingartner 1997; Lazarian & Hoang 2007), i.e., the same direction as that by the
DG mechanism. The alignment by RATs can be more eﬃcient if grains have superparamagnetic inclusions (Lazarian & Hoang 2008) or if the pinwheel mechanism is working with RATs
(Hoang & Lazarian 2009).
The eﬃciency of the RATs alignment varies with strength and spectral energy distribution of the radiation ﬁeld, and thus the size of aligned grains should vary accordingly (Draine& Weingartner 1996; Cho & Lazarian 2005). Therefore, this property makes it possible to examine whether the RATs alignment works or not. Observationally, the maximum wavelength
λmax of polarization in or around dark nebulae was shown to be correlated with extinction AV
in the V-band (Whittet et al. 2001; Andersson & Potter 2007). Andersson & Potter (2010)
showed that grain alignment is enhanced by the stellar radiation in the vicinity of a young1

The direction of alignment with RATs will be parallel to the light beam, if the radiation is so strong thatthe rate of precession by the radiative torques is faster than that of the Larmor precession due to internal
magnetic ﬁeld induced by the Barnett eﬀect.

2

81

star HD 97300 in the Chamaeleon I cloud. Whittet et al. (2008) showed that the polarizationeﬃciency pK /τK in the K-band, where τK is optical depth, decreases smoothly with AV beyond
the region where ice mantle feature was detected. This suggests that the alignment eﬃciency
is not directly related to the state of grain surface, as is expected by the ”pinwheel” alignment.
These results agree with the expectations by the RATs alignment. However, it is still not clear
whether RATs alignment works or not in more diﬀuse clouds.
In this paper, we investigate the polarization quantities and the dust temperature forvarious lines of sight. We observed the stars in/around the Pleiades cluster, which is associatedwith diﬀuse reﬂection nebula, with the low-dispersion spectropolarimeter HBS (Kawabata etal. 1999). Dichroic polarization may be enhanced around those objects, because intense stellar
radiation may make grains align better if the RATs alignment works. The contents of thispaper are essentially the same as those in Matsumura et al. (2011), but include a model ofgrains. Manuscripts are also modiﬁed signiﬁcantly.
2.

Observations and Data Reduction

We observed 8 stars in the Pleiades cluster with the low-dispersion spectropolarimeter
HBS (Kawabata et al. 1999) which is attached to the 1.88m telescope in Okayama AstrophysicalObservatory, from October 30 to November 13 in 2008, and from January 13 to 19 in 2009. The
details of observations are presented in (Matsumura et al. 2011). We also use specropolarimetricdata of various stars by Weitenbeck (1999, 2004), which are similar to our HBS data. The dataof HD 29647 (Whittet et al. 2001) and HD 38087 (Serkowski et al. 1975) are used in addition
to those from Weitenbeck (1999). Also used are the data for high latitude clouds MBM 30 and
MBM 20 (LDN 1642) by Seki & Matsumura (1996).
Extinction Aλ and ratio of total to selective extinction RV for each star are cited from
Fitzpatrick & Massa (2007) when available from their 328 extinction curves. Otherwise, Aλ is
estimated with Fitzpatrick (2004) who showed typical extinction curves for RV = 2.1 − 5.5. For
the stars in the Pleiades cluster2 , we adopt RV = 3.2, the same value as for HD 23512, and thecolor excess is cited from Crawford & Perry (1976). Since RV values of HD 14889 in MBM 30
and SAO149760 in MBM 20 are not available, we set RV = 2.7 or 3.2 as typical values in highlatitude clouds (Larson & Whittet 2005).
Schlegel et al. (1998) deduced dust temperature and other quantities with the COBE
and IRAS data, on the assumption of λ−2 emissivity in the infrared. Their data is homogeneous
all over the sky, so we use their temperature Tdust to examine the radiative environment aroundthe grains in each line of sight.
2

We derived RV values for each star in Matsumura et al. (2011).

3

82

3.

Results

Figure 1 shows the dependence of Aλ , pλ , and pλ /Aλ on wavelength λ for typical examples: HD 23152 and HD 24368 in the Pleiades cluster, HD 29647 in the Taurus cloud, and HD
38087 associated with reﬂection nebula VDB 57. The values of log(pλ /Aλ ) decrease linearlywith inverse wavelength 1/λ in the optical region (Figure 1c), because both log(Aλ ) and log(pλ )are convex against 1/λ (Figures 1a and 1b). We thus make linear ﬁtting in λ = 0.4 − 0.8µm
with the equation:ln(pλ /Aλ ) = ln α − β(1/λ − 1/0.55µm),(1)< 3%/mag forwhere α and β are parameters. The parameter α corresponds to pV /AV , and is ∼interstellar polarization (Serkowski et al. 1975).
Figure 2 presents relations between RV and λmax . Here λmax is the wavelength of maximum polarization pmax , and is obtained with the expression by Serkowski et al. (1975):pλ = pmax exp(−K ln2 (λ/λmax )),(2)

where K = 1.66(λmax /1µm) + 0.01 (Whittet et al. 1992). Long dashed line in Figure 2 showsthe empirical relation between RV and λmax by Serkowski et al. (1975). Most of the stars in
Figure 2 are present near this line except for HD 29647 and HD 38087.
We draw contours of β=0.4, 0.8, and 1.2 µm in Figure 2 on the assumption that polarization pλ follows equation (2), and that extinction Aλ varies as in Fitzpatrick (2004) in
λ = 0.4 − 0.8µm. The contour of β = 0.8µm in Figure 2 is almost coincident with the empirical< 3.5.
line by Serkowski et al. (1975) in RV ∼4.

Discussion

4.1.

Grain Model: Interpretation of βIn this model3 , we assume oblate grains (axial ratio is 2:1) composed of astronomicalsilicate (refractive index is 1.7+0.03i) as a ﬁrst step, although the spheroidal shape is very symmetrical and the radiative torques do not work provided that the incident light is unpolarized.
More complicated shapes, e.g., model in Wurm & Schnaiter (2002), should be necessary, but it
is beyond the scope of this paper. Here we perform light scattering calculation with Fredholm
integral equation method (Matsumura & Seki 1991; Matsumura & Seki 1996; Matsumura &Bastien 2003; Bastien & Matsumura 2005; Matsumura & Bastien 2009), and evaluate extinction and polarization to study the property of β. Following Kim et al. (1994), we assume thesize distribution n(a), where a is a radius of sphere with the same volume, as
3

The modelling is not included in Matsumura et al. (2011). We deleted it according to an anonymous referee’scomments. The parameter β in this paper is diﬀerent from the one often used to characterise the spectrumof thermal emission from dust in far IR. Please do not be confused!

4

83

Aλ [mag.](a) Aλ100

RV=3.2

HD380875

10–1

HD29647 (x 0.2)5pλ [%]

52

1(b) pλHD38087 (x1.5)

HD23512
HD29647 (x0.7)

0.5

HD24368

pλ/Aλ [%/ mag.]

10
5HD24368(x1.5)

HD38087 (x3)

2

HD23
512
(x0.7
)

1
0.5
0.2

(c) pλ/AλHD29647
1

1/λ [µm–1]

2

3Fig. 1. Wavelength dependence of (a) extinction Aλ , (b) fractional polarization pλ , and (c) polarizationeﬃciency pλ /Aλ . Extinction curves in (a) are for HD 29647 (short dashed line) and HD 38087 (long dashedline) from Fitzpatrick & Massa (2007), while that for HD 23512 and HD 24368 (solid line) is the curve of
RV =3.2 from Fitzpatrick (2004). In (b) and (c), asterisks show HD 38087 (Serkowski et al. 1975), ﬁlled
circles HD 23512 (HBS), small open circles HD 24368 (HBS), and large open circles HD 29647 (Whittet
et al. 2001). Long dashed lines in (b) and (c) show HD 29647 and HD 38087 by Weitenbeck (1999).

5

84

λmax [µm]0.80.18µmamin2=0.12µm
µm
β=0.8

29647

µm
1.2
=
β
0.6

38087
210121
0.4

µm
0.4
=
β

amin2=0.06µm

4
RV

6

2

Fig. 2. Relations between RV and λmax . Filled circles show two stars in the Pleiades cluster, i.e., HD
23152 (upper) and HD 24368 (lower), open circles Weitenbeck (1999, 2004), triangles Seki & Matsumura(1996), squares Whittet et al. (2001), and asterisk shows HD 38087 by Serkowski et al. (1975). Long dashedline is the empirical relation RV = 5.6(λmax /1µm) by Serkowski et al. (1975). Solid lines are contours for
β=0.4, 0.8, and 1.2µm. Dashed chains and short dashed line are results of model calculations (see text).

n(a) ∝ aγ exp(−a/amax ),(3)

where γ is a power and amax correspond to a parameter for upper cutoﬀ. For lower cutoﬀ, weset as amin = 0.005µm. The size distribution of aligned grains n(a)align is assumed to be diﬀerent
from n(a), and the expression by Mathis (1986) is used:n(a)align ∝ n(a)(1 − exp(−(a/amin2 )3 )),(4)

where amin2 is a parameter that deﬁnes the typical smallest size of aligned grains. Note that theparameters amax and amin2 are not sharp cutoﬀs for size distributions, but appear in the above
equations, while the parameter amin is a ”sharp” cutoﬀ. Although the size distribution assumed
here is simple, our model shows relations between polarization quantities and parameters forsize distribution, and is helpful to interpret the observational results.
We can obtain the average interstellar values of RV = 3.1 and λmax = 0.53µm, if weset (amin2 , amax , γ) = (0.12µm, 0.12µm, −2.75)4 . We draw one of the dashed chains in Figure
2, labeled as amin2 = 0.12µm, with changing only γ from −3.5 to −2, but (amin2 , amax ) are
kept as (0.12µm, 0.12µm). When only amax varies from 0.09 to 0.24µm, while (amin2 , γ) =(0.12µm, −2.75), a short dashed line is drawn. Those lines are almost coincident with the lineof β = 0.8, showing that β is not sensitive to the variation of γ or amax . However, if we change
amin2 , signiﬁcant variation in β is found: β is ∼ 0.4µm for amin2 = 0.06µm, while it is ∼ 1.2µm
for amin2 = 0.18µm (Figure 2). For our model, we can writeβ ≃ 6.6 × amin2 .
4(5)

It is just a coincidence that amin2 = amax , and this does not show any physical meanings.

6

85

β [µm]

1.2

29647

1

192001 23512 193322

0.8

27Tau

0.6

210121

38087

0.4
14

16

18
20
Tdust [K]22

Fig. 3. Correlation between Tdust and β. Symbols are same as those in Figure 2, and upper small ﬁlled
circle shows 27 Tau, and lower one 19 Tau. Double circle is the average of interstellar medium. Resultof linear ﬁtting for 12 stars is drawn as dashed line. Upper dashed line presents model prediction foramin2 = 5alower , and lower one for amin2 = 3alower .

The dependence of β on the parameter amin2 can be interpreted as below. Since ln(pλ /Aλ )is linearly correlated with 1/λ, we can express β as
β≃d ln(pλ /Aλ ) ln((pλ1 /pλ2 )/(Aλ1 /Aλ2 ))≃,d(1/λ)
1/λ1 − 1/λ2

(6)

for λ1 ̸= λ2 . If γ or amax varies and large (small) grains are abundant, both pλ and Aλ increasein longer (shorter) wavelength. These variations are almost canceled in equation (6). However,the variation of amin2 has eﬀect on pλ but of course not on Aλ , thus β is much more sensitive
to amin2 than amax .
4.2.

Polarization Eﬃciency and Dust Temperature

We now discuss the correlation between the value of β and grain temperature Tdust
obtained by Schlegel et al. (1998). A weak correlation is found for those 16 stars in Figure 3,
and the correlation coeﬃcient is ∼ −0.3.
Several stars in Figure 3 make the correlation worse, i.e. 27 Tau, HD 23512, HD 192001,
and HD 193322. If we exclude those stars, the correlation coeﬃcient is ∼ −0.8 for the rest of
12 stars. The temperature for those outliers may be overestimated, because hotter regions arepresent near or towards them. The line of sight to HD 23512 passes a small molecular cloudthat is present in the intracluster of Pleiades (Breger,1986). The regions near O-type stars insky show high Tdust (Schlegel et al. 1998), and this is the case for HD 192001 and HD 193322.
Using their equation (5) in Cho & Lazarian (2005) on RATs theory, with typical values
for physical quantities in interstellar space, we can express the smallest size alower of aligned
grains as7

86(alower /1µm) = 2.8 × 10−2 × (Tdust /18K)−2 ,(7)

for dust temperature Tdust . The derived value of alower ≃ 0.028µm for Tdust = 18K is smaller
than the value amin2 ∼ 0.12µm expected from observation (Section 4.1). Nevertheless, if we
use equations (5) and (7), we can draw in Figure 3 upper and lower short dashed lines foramin2 = 5alower and amin2 = 3alower , respectively. Those lines follow the observations well, i.e.,the variation of β against temperature can be explained with the RATs alignment theory.
Since the data in Figure 3 contain regions of various temperature, i.e., the Taurus dark cloud,
reﬂection nebulae, etc., our results suggest that the alignment by RATs is ubiquitous in theinterstellar space.
This work was supported by the Thesis Supporting Program at Okayama AstrophysicalObservatory of NAOJ, NINS (No.08A-S01, P.I. was Youko Kameura), and by the Kagawa
University Specially Promoted Research Fund (FY2008). We are grateful to the staﬀ members
at Okayama Astrophysical Observatory and the HBS team for their support during the observations. This work has made use of the SIMBAD database, operated at CDS, Strasbourg,
France.
References
Andersson, B.-G., & Potter, S. B. 2007, ApJ, 665, 369
Andersson, B.-G., & Potter, S. B. 2010, ApJ, 720, 1045
Bastien, P., & Matsumura, M. 2005, Astronomical Polarimetry: Current Status and Future Directions,
343, 149
Breger, M. 1986, ApJ, 309, 311
Cho, J., & Lazarian, A. 2005, ApJ, 631, 361
Crawford, D. L., & Perry, C. L. 1976, AJ, 81, 419
Davis, L., Jr., & Greenstein, J. L. 1951, ApJ, 114, 206
Dolginov, A. Z., & Mitrofanov, I. G. 1976, Ap&SS, 43, 291
Draine, B. T., & Weingartner, J. C. 1996, ApJ, 470, 551
Draine, B. T., & Weingartner, J. C. 1997, ApJ, 480, 633
Fitzpatrick, E. L. 2004, Astrophysics of Dust, 309, 33
Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320
Hoang, T., & Lazarian, A. 2009, ApJ, 695, 1457
Jones, R. V., & Spitzer, L., Jr. 1967, ApJ, 147, 943
Kawabata, K. S., et al. 1999, PASP, 111, 898
Kim, S.-H., Martin, P. G., & Hendry, P. D. 1994, ApJ, 422, 164
Larson, K. A., & Whittet, D. C. B. 2005, ApJ, 623, 897
Lazarian, A. 2007, J. Quant. Spectrosc. Radiat. Transfer, 106, 225
Lazarian, A., & Draine, B. T. 1999, ApJL, 520, L67
Lazarian, A., & Draine, B. T. 1999, ApJL, 516, L37

8

87

Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910
Lazarian, A., & Hoang, T. 2008, ApJL, 676, L25
Markkanen, T. 1977, A&A, 56, 83
Mathis, J. S. 1986, ApJ, 308, 281
Matsumura, M., & Bastien, P. 2009, ApJ, 697, 807
Matsumura, M., & Seki, M. 1991, Ap&SS, 176, 283
Matsumura, M., & Seki, M. 1996, ApJ, 456, 557
Matsumura, M., Seki, M., & Kawabata, K. 1999, AJ, 117, 429
Matsumura, M., & Bastien, P. 2003, Grain Formation Workshop 2003, 23, 7
Matsumura, M., Kameura, Y., Kawabata, K. S., et al. 2011, PASJ, 63, L43
Purcell, E. M. 1979, ApJ, 231, 404
Roberge, W. G., & Lazarian, A. 1999, MNRAS, 305, 615
Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
Seki, M., & Matsumura, M. 1996, Polarimetry of the Interstellar Medium, 97, 168
Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, ApJ, 196, 261
Weitenbeck, A. J. 1999, Acta Astronomica, 49, 59
Weitenbeck, A. J. 2004, Acta Astronomica, 54, 87
Whittet, D. C. B., Gerakines, P. A., Hough, J. H., & Shenoy, S. S. 2001, ApJ, 547, 872
Whittet, D. C. B., Martin, P. G., Hough, J. H., Rouse, M. F., Bailey, J. A., & Axon, D. J. 1992, ApJ,
386, 562
Whittet, D. C. B., Hough, J. H., Lazarian, A., & Hoang, T. 2008, ApJ, 674, 304
Wurm, G., & Schnaiter, M. 2002, ApJ, 567, 370

988Laboratory analogy of amorphous enstatite fine grain formationand crystallization
C. Kaito, Y. Saito, C. Koike, T. Miyashita, T. Fujisaki
Department of Physics, Ritsumeikan University, Kusatsu-shi, Shiga 525-8577, Japan

Abstract
Amorphous enstatite (MgSiO3) grains were produced by the simultaneous
evaporation of Mg and SiO vapor by up and down double heater method in Ar gas
pressure of 10 Torr. Produced particles were mainly mixture of MgO crystallites andamorphous MgSiO3 structure. High resolution electron microscopy showed the
crystallites of MgO and MgSiO3 crystallites less than 10 nm. Crystallization
temperature of MgSiO3 amorphous structure haｓ been determined to be 600℃ by
direct observation by HRTEM in vacuum.

Introduction
In a previous paper,1) we reported the laboratory analogy of crystallize enstatite
grain formation in RF plasma field. The produced MgSiO3 grain size was less than 50
nm and covered with amorphous β-cristobalite phase of SiO2. The growth of Mg
particle on the SiO particle growth stream in RF plasma field was characteristic
difference than the formation of Mg2SiO4 phase particle by the use of parallel double
heaters.2) Submicron-sized silicate grains are present in the circumstellar outflow aroundoxygen-rich stars and within the interstellar mechanism. Since the estimated presser ofsolar nebula at 2-3 AU is low,3) experiments on growth mechanism in vacuum areimportant for composition with condensation, evaporation, melting and crystallization
processes observed in the primitive solar nebula.4,5) A new method to produce thecompound by the use of reaction between SiO film and Mg ultrafine particle6) showed
the enstatite phase (MgSiO3) by heating at 450-700℃. The growth of the MgSiO3 phase
took place due to the movement of the Mg atoms into the amorphous SiO film.
In the present paper, amorphous enstatite ultrafine grain and the crystallization of
amorphous enstatite ultrafine particles in vacuum of 10-6 Pa have been elucidated.

89

Experimental procedure
The chamber used for sample preparation was a glass cylinder 17 cm in diameterand 33 cm in height. The two evaporation sources were set up and down positions at 2
cm as indicated schematically in Fig.1 (c) in Ar gas at 10 Torr. Typical shape of theparticle formation in Ar gas at 10 Torr is shown in Fig.1 (b). In the present experiment,
two evaporation sources in V-shaped tungsten boats (50 mm length, 2 mm width and 1
mm depth) were set as shown schematically in Fig.1 (c). SiO powder was evaporated
from a tungsten boat heated at 1600℃. The Mg evaporation source was set 2 cm above
the evaporation source of SiO. The Mg powder was also put on a tungsten boat heatedat 600℃. The same shape and evaporation temperature was also indicated in Fig.1 (b).
The evaporation flow due to the heating of the SiO evaporation source wasspontaneously mixed with the Mg smoke indicated in Fig.1 (b). The smoke shapes arealso shown in Fig.1 (a). The collected specimen was examined using a Hitachi-9000
electron microscope.

Fig.1 (a) smoke produced at 10 Torr Ar gas by heating 1600℃. (b) Typical increment of thetemperature at various height by heating SiO at 1600℃ and Mg at 600℃ in Ar gas at 10 Torr. (c)Smoke shape of SiO (0 mm) and Mg (20 mm) from the tantalum boats. Mg and SiO smokes were
mixture among at 450℃～500℃ temperature region.

90

Crystallization process of the amorphous enstatite was directly observed using a
Hitachi H-9000 HRTEM with a specimen heating holder, which can be heated toapproximately 800℃ in vacuum at 10-6 Pa . Specimen was put on the amorphous holey
carbon film.
Results and Discussion(a) Amorphous MgSiO3 ultrafine particle formation.
If the Mg and SiO powders are evaporated in Ar gas pressure at 10 Torr from theevaporation source of tungsten heater at 600℃ (Mg) and 1600℃ (SiO), the typicalelectron microscopic images of Mg, silicon (Si) and α-crystallites (SiO2) were produced.
Since the gas pressure was 10 Torr, smoke Mg and Si + SiO2 particles were produced.
The Mg crystallites of the order less than 20 nm were oxidized MgO at room
temperature in order to exposure in air. Ultrafine Mg particles were altered as shown in
Fig.2 (a). The diffraction image of MgO shows the micro crystallites on Mg surface.
Fig.2 (b) shows the typical Si + SiO2 with based amorphous diffraction pattern. As
elucidated in a previous paper on direct observation of metamorphism of the siliconoxide grains,7, 8) the produced particles were Si and α-SiO2.
The SiO powder produced by evaporation of SiO powder at the lower heater at
1600℃ was agree with the ultrafine particles heated at the up heater. As indicated in
Fig.1, the Mg and SiO2 particles were mainly coalesce at temperature less than 450℃.
The produced particles were mainly composed of MgO and MgSiO3 amorphous
particles as shown in Fig.3. The amorphous diffraction pattern was corresponding to theintensity rings of MgSiO3. Therefore the produced amorphous particles were mainly
MgSiO3. The dark field image in Fig.3 shows the crystallites of the order of less than 10
nm was produced. As indicated in HRTEM image of amorphous particles in Fig.4, theexistence of crystallites of MgO and MgSiO3 can be seen. Therefore it can be concluded
that MgSiO3 particles are produced by the coalescence growth between Mg and SiO (Si+ SiO2) fine particles.
As elucidated in previous paper on the formation of MgSiO3 crystallites by thereaction between SiO film and Mg particle by heating above 450℃ in vacuum (Saito etal., 1993). The coalescence growth of fine particles of Mg and SiO less than 500℃
becomes imperfect perfect for the MgSiO3 phase. Since the growth speed of Mg and
SiO ultrafine particles with the high speed than the heating in air, the amorphous
MgSiO3 and oxide Mg particles can be produced. Therefore Mg + SiO → Mg + SiO2 +Si → MgO + MgSiO3. As indicated in HRTEM image by the covering of amorphous
SiO2 phase, same particles were covered with amorphous SiO2 grains.

91

Fig.2 Typical Mg (a) and SiO (b) smoke particles at 1 cm and 3 cm from the heater. Mg particle
was oxidized to exposure in air. The image (a) was exposed to the air about 1 min. Evaporated SiO
at 1600℃ is composed of amorphous α-SiO2 and Si cryatallites of a few nm. The coalescencegrowth between. Mg and SiO2 took place as indicated in Fig.3．

Fig.3 Typical electron microscopic image, dark field image and diffraction image produced by
the coalescence growth between SiO and Mg ultrafine particles as indicated in Fig.1. Amorphous
MgSiO3 and MgO fine crystallites with size less than 5nm order crystallites. The small diffraction
pattern shows the small crystallites with 5 nm order MgO and MgSiO3.

92

The IR spectra of the MgSiO3 amorphous grains produced by the present method
have been measured by dispersing the produced particle on KBr crystal as well as theexperiment report on crystalline silicate used the crystalline materials (Chihara et al.,2002). The data was presented in Fig.5. Typical examples on the distribution of MgOand MgSiO3 shown in Fig.4 specimen on KBr crystal are also shown in Fig.5. Nano size
specimen was distributed on the surface of KBr. The data on Oen was appeared as sownin Fig.5. Since MgSiO3 specimen was amorphous and the crystallites of MgSiO3 were

Fig.4 Typical high resolution images from particles in Fig.3. Mg and MgSiO3 crystallites
formation can be seen.

less than 10nm, the IR spectra for a long range were not appeared. In the present
specimen the oxidized MgO crystallites of Mg grain was included, the peak due to the
Mg and O distances in MgO crystallites showed the peak at 14.5 nm as elucidated inprevious paper (Tamura et al. 2003).10) The appearance of forsterite may be due to thelarge SiO2 layer as elucidated in previous paper,１） Mg2SiO4 phase appearance may bemainly due to the process of the sample in KBr.
MgSiO3 and oxide Mg particles can be produced. Therefore Mg + SiO → Mg +SiO2 + Si → MgO + MgSiO3. As indicated in HRTEM image by the covering ofamorphous SiO2 phase, same particles were covered with amorphous SiO2 grains.

93

50nm
Fig.5 IR spectrum collected on KBr crystallites. The produced amorphous particle size of 20 nm
order was distributed on KBr crystal. The formation of MgSiO3 and MgO can be clearly seen. The
existence of Mg2SiO4 may be due to specimen prepared of IR on KBr crystal. MgSiO3 amorphous
particles less than 30 nm were distributed on KbR specimen. KBr specimen was week in theelectron microscope.

The IR spectra of the MgSiO3 amorphous grains produced by the present method havebeen measured by dispersing the produced particle on KBr crystal as well as theexperiment report on crystalline silicate used the crystalline materials (Chihara et al.,2002). The data was presented in Fig.5. Typical examples on the distribution of MgOand MgSiO3 shown in Fig.4 specimen on KBr crystal are also shown in Fig.5. Nano size
specimen was distributed on the surface of KBr. The data on Oen was appeared as sownin Fig.5. Since MgSiO3 specimen was amorphous and the crystallites of MgSiO3 wereless than 10nm, the IR spectra for a long range were not appeared. In the present
specimen the oxidized MgO crystallites of Mg grain was included, the peak due to the
Mg and O distances in MgO crystallites showed the peak at 14.5 nm as elucidated inprevious paper (Tamura et al. 2003).10) The appearance of forsterite may be due to thelarge SiO2 layer as elucidated in previous paper,１） Mg2SiO4 phase appearance may bemainly due to the process of the sample in KBr.
(b) Direct observation on the crystallization of amorphous MgSiO3 particle
The specimen holder on the heating temperature at room temperature to 1073 K
which can be used in previous paper for composition Mg2SiO4 particle (Kamitsuji etal.,) was used for the enstatite ultrafine in direction in Fig.2. Fig.7 shows the typicalelectron microscopic image at 270℃, 400℃ and 600℃ of the same position. Particles
less than 50 nm hardly altered. The small crystal structure difference can be seen
between at 230℃ and 400℃ as well as the low temperature crystallization observed

94

Fig.6 Alteration of grains by heating in air. Crystallites formation took place at 600℃.

Fig.7 High resolution electron microscopic image after heated at 600oC. The image of a particle
indicated clearly shows the growth of amorphous SiO2. The central part of the image clearly
showed the growth of MgSiO3 crystallites. If the SiO layer is larger, the growth of SiO2 layer on
the MgO and MgSiO3 crystallite surface took place.

95

Fig.8 MgSiO3 crystallization more took. If the crystallite particles were heated at 800oC order,
the MgSiO3 crystal growth clearly occurred as indicated diffraction pattern of the particles. The
crystal image growth can be clearly seen in the image of the same position at 600oC, 650oC and
830oC.

Fig.9 The heated specimen at eight hundred order was more crystallized during the coolingprocess. MgSiO3 and MgO crystal growth took placed as seen the diffraction pattern.

Mg2SiO4 amorphous structure (Kaito et al., 2007). The massive appearance of small
crystallites in MgSiO3 grains can be clearly observed by heating at 600℃.The finecrystallites formation less than 10 nm shows the amorphous diffraction. As indicated in
Fig.8, the image at 600℃ is composed of crystallites of the order of 5 nm as seen theimage and dark field image. Some of the particle surface was also covered with SiO2
layer of amorphous structure. The typical examples are shown in Fig.8. The crystalline
MgSiO3 is clearly seen of the particle. As indicated by the dark field image on the

96

produced particles in Fig.7, the crystallites of MgO and MgSiO3 are increased.
Therefore the crystallization took place at 600 ℃ in vacuum of 10-6 Pa. The
crystallization growth by the coalescence took place above 600℃.
Figure 9 shows the result increased 600℃ to 830℃ order, the crystal form werehardly changed. The crystallization took place. The electron diffraction pattern shows
clearly appearance of MgSiO3 crystallites. The electron diffraction pattern of MgO isdiffuse. If the specimen heating system was cut off, the specimen was cooled to the row
temperature at about 2 hours in vacuum system. After to the room temperature, theheated specimen at 800℃ order was more crystallized during the cooling process asshown in Fig.9. If the specimen heated at 600℃ from one hour, the specimen such as
crystallization could not observed. The micro-crystallization took place above 600℃. If
the heating temperature is order of 800℃, the crystallization during the heating stage at
about 600℃ region may be accelerated the crystallization.

Reference
[1]
[2]
[3]
[4]
[5]
[6]
[7]

T. Sato, K. Kamitsuji, M. Shintaku, Y. Kimura, M. Kurumada, O. Kido, H. Suzuki,
Y. Saito, C. Kaito, 2006, Planetary and Space Science 54, 615.
C. Kaito, Y. Ojima, K. Kamitsuji, O. Kido, Y. Kimura, H. Suzuki, T. Sato, T.
Nakada, Y. Saito, C. Kaito, 2003, Mteoritics and Planetary Science 38, 49.
H. Nagahara, I. Kushiro, B. M. Mysen., 1994, Cosmochem. Acta. 58, 1951.
J, A. Nuth, 1996 in The Cosmic Dust Connection, ed. J. M. Greenberg (TheNetherlands : Klumer Academic), 205.
B. O. Mysen, I. Kushiro, 1988, Am. Minerral. 73, 1.
Y. Saito, K. Ohtsuka, T. Watanabe, C. Kaito, 1993, J. Crystal Growth 128, 271.
T. Morioka, S. Kimura, N. Tsuda, C. Kaito, Y. Saito, C. Koike., 1998, Mon. Not. R.

Astron. Soc. 299. 78.
[8] K. Kamitsuji, H. Suzuki, Y. Kumura, T. Sato, Y. Saito, C. Kaito., 2005, Astronomy
and Astrophysics 429, 205.
[9] H. Chihara, C. Koike, A. Tsuchiyama, T. Tachibana, D. Sakamoto, 2002,
Astronomy and Astrophysics, 391, 267.
[10] K. Tamura, Y. Kimura, H. Suzuki, O. Kido, T. Sato, T. Tanigaki, K. Kurumada, Y.
Saito and C. Kaito, 2003, Jpn. J. Appl. Phys. 42, 7489.
[11] C. Kaito, S. Sasaki, Y. Miyazaki, A. Kumamoto, M. Kurumada, K. Yokoyama, M.
Saito and Y. Kimura, 2006, Advances in Geosciences 7, 125.

97

98

99

100

101

Coagulation and Radial Drift of Porous Dust Aggregates in Protoplanetary
Disks
Satoshi Okuzumi1 , Hidekazu Tanaka2 , Hiroshi Kobayashi1 , and Koji Wada3
ABSTRACT
Rapid orbital drift of macroscopic dust particles is one of the major obstacles against
planetesimal formation in protoplanetary disks. We reexamine this problem by considering
porosity evolution of dust aggregates. We apply a porosity model based on recent N-bodysimulations of aggregate collisions, which allows us to study the porosity change upon collision for a wide range of impact energies. As a first step, we neglect collisional fragmentationand instead focus on dust evolution outside the snow line, where the fragmentation has beensuggested to be less significant than inside the snow line because of a high sticking eﬃciency
of icy particles. We show that dust particles can evolve into highly porous aggregates (withinternal densities of much less than 0.1 g cm−3 ) even if collisional compression is taken into
account. Furthermore, we find that the high porosity triggers significant acceleration in collisional growth. This acceleration is a natural consequence of particles’ aerodynamical property
at low Knudsen numbers. Thanks to this rapid growth, the highly porous aggregates are foundto overcome the radial drift barrier at distances less than 10 AU from the central star. Thissuggests that, if collisional fragmentation is truly insignificant, formation of icy planetesimals
is possible via direct collisional growth of submicron-sized icy particles.

Subject headings: dust, extinction — planets and satellites: formation — protoplanetary disks

1. Introduction
Growth of dust particles is a key process in protoplanetary disks. Current theories of planet
formation assume km-sized solid bodies called “planetesimals” to form from dust contained in protoplanetary disks. Being the dominant component of disk opacity, dust also aﬀects the temperatureand observational appearance of the disks.
Theoretically, however, it is poorly understand how the dust particles evolve into planetesimals. One of the most serious obstructs is the radial inward drift of macroscopic aggregates due
1

Department of Physics, Nagoya University, Nagoya, Aichi 464-8602, Japan; okuzumi@nagoya-u.jp

2

Institute of Low Temperature Science, Hokkaido University, Sapporo 060-0819, Japan

3

Planetary Exploration Research Center, Chiba Institute of Technology, Narashino, Chiba 275-0016, Japan

102

–2–
to the gas drag (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977). Because of the gaspressure support in addition to the centrifugal force, protoplanetary disks tend to rotate at subKeplerian velocities. By contrast, dust particles are free from the pressure support, and hence theytend to rotate faster than the gas disk. The resulting head wind acting on the dust particles extracts
their angular momentum and thus causes their drift motion towards the central star. In order to go
beyond this “radial drift barrier,” dust particles must decouple from the gas drag (i.e., grow large)
faster than they drift inward. However, previous work by Brauer et al. (2008) showed that dust
particles finally fall onto the central star unless the initial dust-to-gas mass ratio is considerably
higher than the canonical interstellar value.
This study reexamines this problem by considering a new physical eﬀect: porosity evolutionof dust aggregates. Most previous coagulation models assumed that dust particles grow with afixed internal density, or a fixed porosity (e.g., Brauer et al. 2008). In reality, however, the internal
density of aggregates does change upon collision depending on the impact energy. We simulate
the temporal evolution of the radial size distribution of aggregates using the advection-coagulation
model developed by Brauer et al. (2008). Unlike the previous work, we allow the porosities ofaggregates to change upon collision, depending on their impact energies. To do so, we adopt the“volume-averaging method” proposed by Okuzumi et al. (2009). In this method, aggregates of
equal mass are regarded as having the same volume (or equivalently, the same internal density) ateach orbital distance, and the advection-coagulation equation for the averaged volume is solved
simultaneously with that for the radial size distribution. To determine the porosity change upon
collisional sticking, we use an analytic recipe presented by Suyama et al. (2012) that well reproduces the collision outcomes of recent N-body simulations (Wada et al. 2008; Suyama et al. 2008)as a function of the impact energy. These theoretical tools allow us to study for the first time how
the porosity evolution aﬀects the growth and radial drift of dust aggregates in protoplanetary disks.
Details of our results are presented elsewhere (Okuzumi et al. 2012).

2. ModelWe adopt the minimum-mass solar nebular model of Hayashi (1981) for the disk model. Dust
particles are assumed to be initially compact spheres (monomers) of equal size a0 = 0.1 µm andequal internal density ρ0 = 1.4 g cm−3 , distributed in the disk with a constant dust-to-gas surface
density ratio Σd /Σg = 0.01. Turbulence with the alpha value αD = 10−3 is considered to determine
the vertical profile and collision velocity of aggregates.
We solve the evolution of the radial size distribution of dust aggregates using the method
developed by Brauer et al. (2008). This method assumes the balance between sedimentation andturbulent diﬀusion of aggregates in the vertical direction. We neglect electrostatic and gravitational
interactions between colliding aggregates and assume perfect sticking upon collision. Thus, thecollisional cross section is simply given by σcoll = π(a1 + a2 )2 , where a1 and a2 are the radii of the

103

–3–
colliding aggregates.
The velocity of an aggregate relative to the gas generally depends on the its stopping timet s . To calculate t s as a function of aggregate size, we assume Epstein’s law for a < 9λmfp /4 andStokes’ law for a > 9λmfp /4, where λmfp is the mean free path of gas particles.
The radial drift velocity is given by vr = −2ηvK Ωt s /[1 + (Ωt s )2 ], where 2η is the ratio of
the pressure gradient force to the stellar gravity in the radial direction and vK = rΩ is the Kepler
velocity (e.g., Whipple 1972; Adachi et al. 1976; Weidenschilling 1977). The radial drift speed
has a maximum ηvK , which is realized when Ωt s = 1. In our disk model, η scales with r as
η = 4.0 × 10−3 (r/5 AU)1/2 , and the maximum inward speed ηvK = 54 m s−1 is independent of r.
We also consider the collisional evolution of the porosity of aggregates. To do so, we treat
the volume V = (4π/3)a3 of an aggregate as a variable independent of the mass m. We follow the
temporal evolution of V(r, m) using the volume-averaging method developed by (Okuzumi et al.
2009).
Here we describe our recipe for determining the volume of merged aggregates. In general, thevolume of an aggregate after sticking depends on the impact energy Eimp relative to the energy Erollrequired for one monomer to roll over 90 degrees on the surface of another monomer in contact(Dominik & Tielens 1997; Blum & Wurm 2000; Wada et al. 2007) When Eimp ≪ Eroll , no visible
restructuring occurs upon collision. In this case, the volume after collision is determined independently of Eimp . When Eimp ≳ Eroll , internal restructuring occurs through inelastic rolling among
constituent monomers. In this case, the final volume depends on Eimp as well as on the volumes ofaggregates before collision. Recently, Suyama et al. (2012) have proposed an analytic formula thatcan be used even for arbitrary values of Eimp /Eroll . We use this formula in our simulations. The
rolling energy is assumed to be Eroll = 1.8 × 10−8 erg.
The collision velocity between aggregates is assumed to be driven by Brownian motion, radial
and azimuthal drift, vertical settling, and turbulence. We use the analytic formula of Ormel & Cuzzi(2007) to calculate the turbulence-driven velocity.

3.

Results

To begin with, we show the result of compact aggregation. In this simulation, we fixed theinternal density ρint ≡ m/V of the aggregates to the material density ρ0 = 1.4 g cm−3 , and solved
only the evolutionary equation for the radial size distribution. Figure 1 shows the snapshots ofthe dust surface density per decade of aggregate mass, ∆Σd /∆ log m. At each orbital radius, dust
growth proceeds without significant radial drift until the stopping time of the aggregates reachesΩt s ∼ 0.1 (the dotted lines in Figure 1). However, as the aggregates grow, the radial drift becomes
faster and faster, and further growth becomes limited only along the line Ωt s ∼ 0.1 on the r–m

104

–4–

Fig. 1.— Aggregate size distribution ∆Σd /∆ log m at diﬀerent times t for the compact aggregation model(ρint = 1.4 g cm−3 ) as a function of orbital radius r and aggregate mass m. The dotted lines mark theaggregate size at which Ωt s exceeds 0.1.

plane. Consequently, a significant amount of dust has been lost from the planet-forming regionr ≲ 30 AU within 105 yr.
Now we show how porosity evolution aﬀects dust evolution. Here, we solve the evolution of
aggregate volume simultaneously with the radial size distribution. The result is shown in Figure 2,
which displays the snapshots of the aggregate size distribution ∆Σd /∆ log m and internal density
ρint = m/V (right four panels) at diﬀerent times t as a function of r and m. The left four panels of
Figure 2 show how the radial size distribution evolves in the porous aggregation. At t < 103 yr,
the evolution is qualitatively similar to that in compact aggregation. However, in later stages,
the evolution is significantly diﬀerent. We observe that aggregates in the inner region of the disk(r < 10 AU) undergo rapid growth and eventually overcome the radial drift barrier lying at Ωt s ∼ 1(dashed lines in Figure 2) within t ∼ 104 yr. In the outer region (r > 10 AU), aggregates do drift
inward before they reach Ωt s ∼ 1 as in the compact aggregation model. However, unlike in thecompact aggregation, the inward drift results in the pileup of dust materials in the inner region(r ≈ 4–9 AU) rather than the loss of them from outside the snow line. This occurs because mostof the drifting aggregates are captured before they reach the snow line by aggregates that alreadyovercome the radial drift barrier. As a result of this, the dust-to-gas mass ratio in the inner regions
is enhanced by a factor of several in 105 yr.
The right four panels of Figure 2 show the evolution of the internal density ρint = m/V asa function of r and m. First thing to note is that the dust particles grow into low-density objects
at every location until their internal density reaches ρint ∼ 10−5 –10−3 g cm−3 . In this stage, the

105

–5–

Fig. 2.— Aggregate size distribution ∆Σd /∆ log m (left four panels) and internal density ρint = m/V (rightfour panels) at diﬀerent times t for the porous aggregation model as a function of orbital radius r andaggregate mass m. The dashed lines mark the aggregate size at which Ωt s exceeds unity.

internal density decreases as ρint ≈ (m/m0 )−1/2 ρ0 , meaning that the dust particles grow into fractal
aggregates with the fractal dimension d f ≈ 2 (Okuzumi et al. 2009). At m ∼ 10−4 –10−6 g, thefractal growth stage terminates, followed by the stage where collisional compression becomes
eﬀective (Eimp ≫ Eroll ). In this late stage, the internal density decreases more slowly or is kept at
a constant value depending on the orbital radius.

4.

Condition for Breaking Through the Radial Drift Barrier

Below, we explain why porous aggregates overcome the radial drift barrier in the inner regionof the disk. We do this by comparing the timescale of aggregate growth and radial drift. For
simplicity, we assume that all dust aggregates are similar in size (this assumption is valid when theaggregate size distribution is peaked at m ∼ ⟨m⟩m ). The growth rate of the aggregate mass m at the
midplane is then given by
dm
Σd
= √
A∆v,
dt
2πhd

(1)

where Σd , hd , A, and ∆v are the spatial mass density (at the midplane), scale height, projected area,
and relative velocity of the aggregates. Equation (1) can be rewritten in terms of the grow timescale
tgrow ≡

( d ln m )−1dt√
4 2π hd ρint a=,
3 ∆v Σd

(2)

106

–6–
where we have used that m = (4π/3)ρint a3 and A = πa2 . What we do here is to compare tgrow withthe timescale for the radial drift given by
tdrift ≡d ln r −1
r
=
.
dt
|vr |

(3)

Now we focus on the stage at which the radial drift velocity reaches the maximum value,i.e., Ωt s = 1. At this stage, the dust scale height is given by hd ≈ (2αD /3)1/2 hg , where hg is the√gas scale height (Youdin & Lithwick 2007). In addition, we set ∆v ≈ αD c s since the collisionvelocity between Ωt s = 1 particles is dominated by the turbulence-driven velocity. Substituting
these relations and hg = c s /Ω into Equation (2), we obtain( Σ /Σ )−1 (ρ a)dgint Ωt s =1
tgrow |Ωts =1 ≈ 40
tK ,
(4)
0.01
Σg
where tK = 2π/Ω is the Keplerian orbital period. By contrast, the drift timescale for Ωt s = 1
particles is
(
)−1
1
η
tdrift |Ωts =1 =
≈ 40
tK .
(5)
ηΩ
4 × 10−3
The ratio (tgrow /tdrift )Ωts =1 determines the fate of dust growth at Ωt s = 1. If (tgrow /tdrift )Ωts =1
is very small, dust particles grow beyond Ωt s = 1 without experiencing significant radial drift;
otherwise, dust particles drift inward before they grow. We expect growth without significant drift
to occur if
(t )
1
grow
≲ ,
(6)
tdrift Ωts =1 10
where the threshold value 1/10 has been chosen based on the fact that tgrow is the timescale for mass
doubling while the particles experience the fastest radial drift over one or two orders of magnitude
in mass. Below, we examine in what condition this requirement is satisfied.
The ratio (ρint a)Ωts =1 /Σg depends on the drag regime at Ωt s = 1. For the Epstein drag regime,
it can be shown that (Okuzumi et al. 2012)( Σ /Σ )−1dg
tgrow |Ωts =1 ≈ 30
tK .
(7)
0.01
Hence, the growth condition is not satisfied for the standard disk parameters η ≈ 10−3 and Σd /Σg =0.01, in agreement with the results of our and previous (Brauer et al. 2008) simulations.
The situation diﬀers in the Stokes drag regime. A similar calculation shows (Okuzumi et al.
2012)( Σ /Σ )−1 ( λ )mfpdgtgrow |Ωts =1 ≈ 60
tK .
(8)
0.01
a Ωts =1
Note that the growth timescale is inversely proportional to the aggregate radius, in contrast to that in
the Epstein regime (Equation (7)) being independent of aggregate properties. This means that dust
aggregates break through the radial drift barrier in the “deep” Stokes regime, (a/λmfp )Ωts =1 ≫ 1.

107

–7–

Fig. 3.— Size a (upper panel) and growth timescale tgrow (lower panel) of dust aggregates at Ωt s = 1 as afunction of orbital radius r for internal densities ρint = 1.4 g cm−3 (solid line), 10−2 g cm−3 (dashed line),and 10−5 g cm−3 (dotted line). The minimum-mass solar nebular with the dimensionless diﬀusion coeﬃcient
αD = 10−3 is assumed for the disk model. The thick black line in the upper panel indicates a = 9λmfp /4,
at which the drag law goes from the Epstein regime to the Stokes regime. The thick black line in the lower
panel shows the drift timescale tdrift at Ωt s = 1, which is the minimum value of tdrift and is independent ofρint .

The internal density of aggregates controls the growth timescale through the aggregate size
a at Ωt s = 1. The upper panel of Figure 3 plots a|Ωts =1 for three diﬀerent values of the aggregateinternal density ρint . If dust particles grew into compact spheres (ρint ∼ 1 g cm−3 ), they would
be inside the Epstein regime when Ωt s = 1 in almost entire parts of the snow region (r > 3AU).
However, if dust particles grow into highly porous aggregates with ρint ∼ 10−5 g cm−3 , the particles
growing at r < 30 AU enter the Stokes regime before Ωt s reaches unity. The lower panel of
Figure 3 shows the two timescales tgrow |Ωts =1 and tdrift |Ωts =1 as calculated from Equations (4) and (5),respectively. We see that compact particles with ρint ∼ 1 g cm−3 do not satisfy the growth condition(Equation (6)) outside the snow line, while porous aggregates with ρint ∼ 10−5 g cm−3 do in theregion r ≲ 10 AU. These explain our simulation results presented in Section 3.

108

–8–
5. Summary and Discussion
We have investigated how the porosity evolution of dust aggregates aﬀects their collisional
growth and radial inward drift. We have shown that icy aggregates can become highly porous even
if collisional compression is taken into account. Our model calculation suggests that the internal
density of icy aggregates at 5 AU falls oﬀ to 10−5 g cm−3 by the end of the initial fractal growthstage and then is kept to this level until the aggregates decouple from the gas motion. Stretchingof merged aggregates at oﬀset collisions, which is not taken into account in our porosity model,
could further decrease the internal density (Wada et al. 2007; Paszun & Dominik 2009).
We have also found that the very high porosity triggers significant acceleration in collisional
growth. This acceleration is a natural consequence of particles’ aerodynamical property in the
Stokes regime, i.e., at particle radii larger than the mean free path of the gas molecules. The
porosity (or internal density) of an aggregate determines whether the aggregate reaches the Stokes
regime before the radial drift stalls its growth. Compact aggregates tend to drift inward before
experiencing the rapid growth, while highly porous aggregates are able to experience it over awide range of orbital radii. Furthermore, the rapid growth has been found to enable the aggregatesto overcome the radial drift barrier in inner regions of the disks. Our model calculation shows thatthe breakthrough of the radial drift barrier can occur in “planet-forming” regions, i.e., at distancesless than 10 AU from the central star. The radial drift barrier has been commonly thought to be
one of the most serious obstacles against planetesimal formation. Our result suggests that, if thefragmentation of icy aggregates is truly insignificant (as suggested by Wada et al. 2009), formationof icy planetesimals is possible via direct collisional growth of submicron-sized icy particles even
without an enhancement of the initial dust-to-gas mass ratio.
The porosity evolution of dust aggregates can even influence the evolution of solid bodies
after planetesimal formation. It is commonly believed that the formation of protoplanets begins
with the runaway growth of a small number of planetesimals due to gravitational focusing. The√
runaway growth requires a suﬃciently high gravitational escape velocity vesc = 2Gm/a relative to the collision velocity. Since the escape velocity decreases with decreasing internal density(vesc ∝ m1/3 ρ1/6
int ), it is possible that a high porosity delays the onset of the runaway and thereby affects its outcome. For example, a recent protoplanet growth model including collisional fragmentation/erosion (Kobayashi et al. 2010) suggests that planetesimals need to have grown to > 1021 g
before the runaway growth begins in order to enable the formation of gas giant planets within theframework of the core accretion scenario. Compaction of large and massive aggregates may occur
through static compression due to gas drag or self gravity. To precisely determine when it occurs isbeyond the scope of this work, but it will be thus important to understand later stages of planetarysystem formation.

The authors thank Tilman Birnstiel, Frithjof Brauer, Chris Ormel, Taku Takeuchi, Takayuki
Tanigawa, Fredrik Windmark, and Andras Zsom for fruitful discussions. S.O. also thanks Cornelis

109

–9–
Dullemond and Shu-ichiro Inutsuka for critical comments. S.O. acknowledges support by Grantsin-Aid for JSPS Fellows (22 · 7006) from MEXT of Japan.

REFERENCES
Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756
Blum, J., & Wurm, G. 2000, Icarus, 143, 138
Brauer, F., Dullemond, C. P., & Henning, Th. 2008a, A&A, 480, 859
Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647
Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35
Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836
Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, submitted
Okuzumi, S., Tanaka, H., & Sakagami, M-a. 2009, ApJ, 707, 1247
Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413
Paszun, D., & Dominik, C. 2009, A&A, 507, 1023
Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310
Suyama, T., Wada, K., & Tanaka, H. 2012, ApJ, submitted
Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320
Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296
Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490
Weidenschilling, S. J. 1977, MNRAS, 180, 57
Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (New York: Wiley Interscience), 211
Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588

This preprint was prepared with the AAS LATEX macros v5.2.

110

Three-dimensional structure of Itokawa regolith particles.
Akira Tsuchiyama
Department of Earth and Space Science, Graduate School of Science, Osaka University
1-1 Machikaneyama-cho, Toyonaka, 560-0043, Japan

Particles of Asteroid Itokawa were successfully recovered by the Hayabusa missionof JAXA (Japan Aerospace Exploration Agency) [1-6]. They are not only the firstsamples recovered from an asteroid, but also the second extraterrestrial regolith to havebeen sampled, the first being the Moon, which was sampled by the Apollo and Lunamissions [7]. It is accepted that most meteorites originate from asteroids, as
demonstrated by orbital determination from observed meteorite falls. The materials on
asteroids have been estimated by comparing reflectance spectra between asteroids andmeteorites. Itokawa has an S-type spectrum, which is one of the major types and similar
to those of ordinary chondrites. Remote-sensing observation by the Hayabusa spacecraft[8] indicate that the materials on asteroid Itokawa are similar to specific types ofordinary chondrites (LL5 or LL6 chondrites). However, there is a spectrum discrepancy
between the asteroid and the meteorites, which has been considered to result from space
weathering [9]. However, there is no direct evidence that the surface material of
Itokawa is LL chondrites suffered with space weathering. Itokawa samples allow a
direct validation of the relation between asteroids and meteorites. In addition, theproperties of Itokawa particles allow studies of regolith formation and evolution.
About 2000 recovered particles (<~300 µm and most of them are <~10 µm) havebeen identified as having an Itokawa origin by curation at JAXA. About fifty particles
were allocated for preliminary examination (PE) [1] and their results were published[1-6,10,11]. We examined three-dimensional (3D) structures of Itokawa particles using
non-destructive x-ray micro-tomography as the first analysis in a sequential PE
analytical flow [2]. The purposes are to understand their materials in comparison withmeteorites and the 3D shape features in connection with regolith formation andevolution. The examination has additional technical purpose to provide information forthe design of later destructive analyses. This is one of the key features of the Hayabusa
PE strategy.

111

Forty particles (~30 to 180 µm) were imaged at BL47XU of SPring-8 with effective
spatial resolution of ~200 or ~500 nm [2] using absorption imaging tomography [12,13].
Imaging at two X-ray energies of 7 and 8 keV made identification of minerals in CT
images possible (Fig. 1) since the K-absorption edge of Fe (7.11 keV) is present
between the two energies (Fig. 2) (analytical dual-energy micro-tomography [14]). The
Fe/Mg ratios of olivine, low-Ca pyroxene and high-Ca pyroxene, and the Na/Ca ratio of
plagioclase can be also obtained using this method. A successive set of 3D CT images,
which shows quantitative 3D mineral distribution, was obtained for each particle.

Figure 1. CT images of an Itokawa particle (RA-QD02-0024). (A) 7 keV. (B) 8 keV.
Minerals can be identified by comparing their brightness at two energies, whichcorresponds to linear attenuation coefficients of minerals. Ol: olivine, LPx: low-Ca
pyroxene, HPx: high-Ca pyroxene, and Pl: plagioclase.

112

Figure 2. A two-dimensional histogram of linear attenuation coefficient (LAC) values at
7 and 8 keV for an Itokawa particle (RA-QD02-0024) (original figure in [3] waschanged). Solid symbols correspond to the mean chemical compositions of the minerals
in Hayabusa samples [1]. Fo: forsterite, En: enstatite, Di: diopside, Ab: albite, An:anorthite. The other abbreviations are the same as those in Figure 1. Numbers alongolivine, pyroxene and plagioclase solid solutions are the forsterite, enstatite,wollastonite and albite contents (in mol.%).
The tomography gave the total volume of 4.2x106 µm3, which corresponds to asphere of ~200 µm in diameter (total mass of ~15 mg). The mineral assemblage andtheir abundances of the whole samples are similar to those of LL chondrites (Fig. 3).
Slightly larger amount of olivine and smaller amounts of high-Ca pyroxene, troilite and(Fe,Ni) metal of the Itokawa particles can be regarded as statistical errors due to picking
up the small sample amount [15]. Other PE data, such as the elemental and oxygenisotopic compositions of minerals [1,3,4], also indicate LL chondrites.

Figure 3. Modal mineral abundances (in vol.%) of the whole Itokawa particles and theaverages of LL4-6, L4-6 and H4-6 chondrites. Tr: troilite. The other abbreviations are
the same as those in Figure 1.
3D structures of the most particles have highly-equilibrated textures by thermal
metamorphism (Figs. 4A-C), which correspond to the petrologic type of 5 and/or 6 inordinary chondrites (LL5 and/or LL6), while some of them (~10%) haveless-equilibrated textures (LL4) (Fig. 4D). SEM observation and compositional
homogeneity of minerals show the same results [1]. Fe-rich nanoparticles observed inthin surface layers (<60 nm) of Itokawa particles show the evidence of space

113

weathering [5]. These results showed that the Itokawa surface material is consistent
with LL chondrites suffered by space weathering as expected and brought an end to theorigin of meteorites. In future detailed studies on 3D micro-textures of LL chondrites,
including voids and micro-cracks, using the same resolution as the micro-tomography
used in the PE of the Itokawa particles, should be made systematically to compare
Itokawa samples with LL chondrites because Itokawa samples are so small that the
textures cannot be compared with those of LL chondrites using conventional
observations.

Figure 4. Representative CT images of Itokawa particles [2]. (A) Sample
RA-QD02-0063. (B) RA-QD02-0014. (C) RA-QD02-0042. (D) RA-QD02-0048.
Concentric structure is a ring artifact. Bright edges of particles and voids are artifacts
resulting from refraction contrast. CP: Ca phosphate and Meso: mesostasis. The otherabbreviations are the same as those in Figures 1 and 3.
External shapes of the Itokawa particles were extracted from the 3D CT images.
Sphere-equivalent diameters and three-axial lengths of each particle were obtained fromthe shapes using image analysis [2]. Cumulative size distribution of the particles has thelog-slope of ~-2. This is consistent with dominance of ~cm size regolith in the Itokawasmooth terrain, which was proposed by remote-sensing observation by the Hayabusa

114

spacecraft [16], by comparing with the size distribution of boulders (0.1-5 m) of thelog-slope of ~-3 [17]. The 3D shape (three axial length ratios) distribution cannot be
distinguished from that of fragments produced by impact experiments [18] (Fig. 5A),
showing that the Itokawa particles are consistent with impact fragments. No particles
showing melting were observed, indicating relatively low-impact velocities similar totypical relative impact velocities among asteroids (~5 km/s, in contrast to lunar regolith,
where agglutinates formed by high-impact velocities of >~10 km/s [7]). Most of theparticles have sharp edges (Fig. 6A), while some (~25%) have rounded edges (Fig. 6B).
They were probably formed from particles that were originally more angular byabrasion as grains migrate during impacts (Fig. 7). The spherical shapes of lunar
regolith particles (Fig. 5B) may be due to their longer residence time (~1 billion years)
in the regolith than Itokawa regolith (<~3 million years [6]). Systematic studies on
regolith shapes together with micro- and nano-structures of particle surfaces using
FE/SEM and TEM and solar wind noble gas analysis will give comprehensive
understanding of processes, which should be called “space weathering” in a broad sense,on celestial bodies without atmosphere not only of small steroids, such as Itokawa, butalso of the Moon.

Figure 5. The 3D shape distributions of (A) Itokawa particles and (B) lunar regolith [2].
Fragments of impact experiments are also shown. Large circles in (A) show particles
with rounded edges.

115

Figure 6. Stereograms showing the 3D external shapes of Itokawa particles of (A)RA-QD02-0023 (232 × 232 × 203 µm) and (B) RA-QD02-0042 (112 × 112 × 93 µm)[2].

Figure 7. A schematic illustration of possible abrasion process of regolith particles on
Itokawa by granular processes induced by seismic vibration due to impact of cm-size
meteors.

116

Acknowledgement
The author is grateful to team members of the tomography, the PE, and the
Hayabusa project.
References[1] Nakamura, T. et al., 2011. Science 333, 1113-1116.
[2] Tsuchiyama, A. et al., 2011. Science 333, 1125-1128.
[3] Yurimoto H. et al., 2011. Science 333, 1116.
[4] Ebihara M. et al., 2011. Science 333, 1119.

[5] Noguchi T. et al., 2011. Science 333, 1121.
[6] Nagao K. et al., 2011. Science 333, 1128.

[7] G. H. Heiken et al., 1991. Eds., Lunar Sourcebook (Cambridge Univ. Press,
Cambridge).
[8] Abe M. et al., 2006. Science 312, 1334.
[9] B. Hapke, 2001, J. Geophys. Res. 106, 10039.
[10] Naraoka H. et al., 2011. Geochem. J. in press.

[11] Kitajima F. et al., 2011. Lunar Planet. Sci. XLII, 1855. pdf. 28.
[12] K. Uesugi et al., 2006. Proc. SPIE 6318, 63181F.
[13] A. Takeuchi et al., 2009. J. Phys. Conf. Ser. 186, 012020.
[14] Tsuchiyama A. et al., submitted to Geochim. Cosmochim. Acta.
[15] Nagano, T. 2012. Abstract in Lunar and Planetary Science Conference, IVIII.
[16] H. Yano et al., 2006, Science 312, 1350.
[17] T. Michikami et al., 2008, Earth Planets Space 60, 13.
[18] F. Capaccioni et al., 1984, Nature 308, 832.

117

An synthesis experiment of GEMS analogue material produced by thermal plasma.
Junya Matsuno1, Masahito Yagishita2, Takayuki Watanabe2, Akira Tsuchiyama1
1)

Department of Earth and Space Science, Graduate School of Science, Osaka University
1-1 Machikaneyama, Toyonaka, Osaka 556-0043, Japan2)

Department of Environmental Chemistry and Engineering, Tokyo Institute of Technology
4259 G1-22, Nagatsuta, Yokohama, Kanagawa 226-8502, Japan

1. Introduction
Anhydrous IDPs (interplanetary dust particles), which is considered to have cometary origin, have alarge amount of amorphous materials called GEMS (glass with embedded metal and sulfides). They arespherical materials of SiO2-rich silicate glass typically a few 100 nm in diameter that include small(typically 10−50 nm in size) and rounded grains of Fe, Ni metal and sulfides. There are two origins
proposed for GEMS: (1) an interstellar origin (Bradley & Dai 2004) and (2) the solar system. (2) they
were formed in the solar nebula as non-equilibrium condensates from residual Si-rich gas after
condensation of MgO-silicate crystals (Keller & Messenger 2011). It also supports their solar system
origin that 94–99% of GEMS grains have oxygen isotopic compositions that are indistinguishable fromterrestrial materials and carbonaceous chondrites.
Their origin and formation process is still controversial although there are a number of studies on
GEMS in IDPs based on astronomical observations using IR spectrum (Bradley et al. 1999) andlaboratory analyses, such as, microstructures, isotopic and elemental compositions) (e.g., Bradley & Dai
2004, Keller & Messenger 2011). It is important to reproduce GEMS-like materials in a laboratory tounderstand conditions for GEMS formation process. Especially, we focused on whether or not amorphoussilicate spheres that contain nano-grains of metal can be formed by condensation from a Si-rich gas in thisstudy.
2. Experimental method
Induction thermal plasma (ITP) method (Fig. 1) was used in the present experiments. ITP canprovide ultra-high temperatures (~ 10,000 K) to evaporate a starting material immediately. Then, the gas
is quenched rapidly with the cooling time scale of 104 – 105 K/s to form nanoparticles, which are usuallyin an amorphous state. Powders of regent-grade MgO, metallic Fe and SiO2 were mixed together with the
GEMS mean composition (Mg/Si = 0.65, Fe/Si = 0.56) for preparing the starting material, which is ananalogue to Si-rich gas in the early solar system. Mg, Si, Fe and O were taken into consideration as major
elements in solid materials in the solar system for simplicity. S was excluded because of the experimental
difficulty.
The ITP experiments were made under an Ar-He atmosphere at 1 bar at Watanabe laboratory ofTokyo Institute of Technology (TP12010, JEOL). The mixed Ar and He gas is used as the sheath gas andthe flow rate is 60 and 5 L/min, respectively. 6.5 L/min of Ar gas was used for carrier gas and 5 L/min of
Ar was used for inner gas. The plasma power was 30 kW and current frequency was 4 MHz. It should be
noted that the oxygen source came only from the starting material of MgO and SiO2 (Yagishita et al.
2011).
Run products attached on the inner chamber walls of the furnace were collected and they were
exposed to CuKα radiation (λ=1.54056 Å) using a powder X-ray diffraction system (XRD; Ultimate,
Rigaku) at Tokyo Institute of Technology. High spatial resolution images were obtained by using atransmit electron microscope (TEM: JEM-2001, JEOL; JEM-2010F, JEOL) at Tokyo Institute ofTechnology with the accelerating voltage of 200 keV. TEM tomography was also applied to one sample(Fig. 3a). Some run products were embedded in an epoxy resin, and the cross section was observed using
a field emission scanning electron microscope (FE-SEM: JSM-7001F, JEOL) at Osaka University. The
chemical compositions were measured using energy-dispersive X-ray spectroscopy (EDS) equipped with
the TEM (NORAN Instrument) or the FE-SEM (Inca, Oxford).

118

3. Results
The powder XRD pattern of the run product (Fig. 2) shows two sharp peaks at ~45˚ and ~65˚ with abroad halo peak around 27˚ in 2θ. The former peaks were identified as α-iron and the halo as anamorphous silicate. No crystalline silicates, such as forsterite, pyroxene and silica mineral, have been
detected. Thus, amorphous silicate was formed directly from high-temperature gas by very rapid
condensation and cooling. Bright field images of TEM (Figs. 3ab) show that the run product wascomposed of numerous spherical grains (typically ~50 nm in diameter). The TEM tomography suggeststhat each grain has an iron core (~20 nm in diameter) embedded in an amorphous silicate. Figure 4 showssize distributions of amorphous silicate grains of ≥ 30 nm in a diameter and metal grains of ≥ 10 nm in adiameter, respectively. Both size distributions are very broad and no clear peaks were seen. The bulk
chemical composition of the run product was almost consistent with that of the starting material although
the amount of MgO was slightly smaller than that of the starting material (Fig. 5).
4. Discussions
4.1 Formation process of the run product
The texture of the run product, where each metallic iron grain seems to be included always in each
amorphous silicate sphere, strongly suggests that the iron grain condensed first from high-temperature gas
by homogeneous nucleation, and a silicate melt (or the amorphous silicate) nucleates heterogeneously on
the metal grin and covered the iron grain by growth. If the silicate melt condensed first, metallic ironshould nucleate on the silicate melt, and then the silicate melt cover the metal.
The possibility of the heterogeneous nucleation can be checked by considering the surface tensions ofmetallic iron and amorphous silicate as shown in Figure 6. The surface tensions of iron, γ i, and silicate
melt with chondritic composition, γ c, are estimated as 1.75 and 0.4 J/m2, respectively (Table 1) (Dyson
1963; King 1951). However, the surface tension of a boundary between chondrule melt and iron, γ b, was
not known. Uesugi et al. (2008) estimated the gb value of ~1.9 J/m2 from chondrule textures in thinsections, where a small iron globule is attached onto the chondrule surface, by assuming the surface
tension equilibrium. From these values, the total surface energies for the cases where the iron is totallyembedded in the silicate, Et,IN, and the iron and the silicate separated each other, Et,OUT (Fig. 6 and Table1) by using the following equation,
𝐸!,!"#   = 4𝜋[𝛾! 𝑟! ! +    𝛾! 𝑟! ! −    𝑟! ! !/! ](3.1)

𝐸!,!"   = 4𝜋(𝛾! 𝑟! ! +    𝛾! 𝑟! ! )(3.2)

where, gs is the surface energy of the silicate melt (or amorphous silicate),ri is the radius of the metal grain, and rs is the radius of the silicate grain. The values of Et,OUT (2.1×10-14 J)
is a little smaller than the value of Et,IN (2.2×10-14 J), if we take rs = 50 nm and ri = 20 nm and assume
gs=gc (Table 1). Et,IN > Et,OUT indicates that the heterogeneous nucleation of the silicate melt on the iron (orof the metal on the silicate melt) and covering the silicate melt on the iron are not possible. Strictly
speaking, the values of surface tensions depend on chemical compositions and temperatures. Therefore
this estimation should be reconsidered using the surface tensions of silicate melt and boundary of silicate
melt and iron corresponding to the present chemical composition.
Another clue to understand the formation process of the run product is experiments using different
metals other than iron in the same silicate. Different textures might form due to different condensation
temperatures and different surface tensions of the metals. Additional experiments of the metal – silicate
system is on going (Yagishita, private communication).
4.2 Grain formation condition: comparison with ITP and astronomical environments
It is most likely that GEMS grains formed by condensation from a high-temperature gas either incircumstellar regions of evolved stars or the protoplanetary disk of the solar system. Yamamoto &Hasegawa (1977) theoretically formulated homogeneous nucleation and growth of dust grains from a gas,and proposed a non-dimensional parameter for the condensation, Λ, by considering the concentration in agas, c1(0), and cooling time scale of the gas, tT (Table 2). The values of Λ are given for someastronomical environments, such as presolar nebula at 0.1 AU (~3×109) or around AGB stars (0.9~90)

119(Table 2).
Parameters in the present ITP experiments are listed in Table 3. The value of τT was estimated to be~1×10-1 sec based on a temperature profile in the ITP furnace and a flowing gas velocity obtained by anumerical calculation (Watanabe, private communication). Then, Λ in the present experiments wasestimated to be ~4×103 (Table 3). This value is not the same as but not extremely different from those in
the astronomical environments.
Two dashed curves in Figure 4 show size distributions of the silicates and iron calculated by thecondensation theory of Yamamoto & Hasegawa (1977) with the parameters in Table 3 andnon-dimensional parameter for surface tension, μ (silicate ~20 and iron ~10). The calculated size
distributions have sharp peaks at ~42 nm and ~65 nm for the silicate and the metal, respectively. The
peak position of the silicate is roughly similar to that in the experiments, but that of the iron is large. As
this grain formation theory was based on homogeneous nucleation, the calculated size distribution of theiron grain may indicate that the iron in run products did not form by homogeneous nucleation.
4.3 Application to GEMS formation
The textures of the run product are similar to that of GEMS although GEMS is composed of morethan one metal grain. The GEMS texture can form by sintering a number of silicate spheres including a
metal grain. Sintering experiments should be made and estimate temperature for the sintering. It should
be noted that the present experiment does not contain sulfur although FeS grains are included as well as
metal grains in GEMS. FeS grains are distributed in GEMS’s rim, and FeS may form by reaction of ironwith surrounding S-bearing gas in the primordial solar nebula (Keller & Messenger 2011).
Solar system origin of GEMS is proposed based on rare oxygen isotopic large anomalies of GEMS
(Keller & Messenger 2011). However, if GEMS is a mixture of primary grains of a few tens nm in size,
exchange of oxygen atoms between the primary grains and a surrounding gas, which contains large
amounts of H2O and CO molecules, might occur even at low temperatures, and the oxygen isotopic
anomalies disappeared in most of GEMS. Therefore, the rare oxygen isotopic large anomalies of GEMS
might not be evidence of the solar system origin. Primary grains of GEMS might originally form aroundevolved stars by condensation, was transferred to interstellar region, and irradiate by cosmic ray. Then,the primary grains were incorporated into the primordial solar nebula, suffered by oxygen isotope
exchange with surrounding gas, and accumulated into GEMS, and finally some of iron grains near the
GEMS surfaces were sulfidized. More experimental and theoretical work, such as sintering experiments
and estimation of oxygen isotope exchange rate, are required to confirm the above GEMS formation
model.
Acknowledgement
We thank Dr. D. Lu and Mr. N. Hatakeyama for the analysis by the transmission electron microscopy.
References
Bradley, J. P., et al. 1999, Science, 285, 1716
Bradley, J. P., & Dai, Z. 2004, ApJ, 617, 650
Dyson, B. F. 1963, Transact. Metallurgical Soc. of AIME, 277, 1098
Keller, L. P., & Messenger, S. 2011, Geochim. Cosmochim. Acta, 75, 5336
King, T. B. 1951, J. Soc. Glass Tech., 35, 241
Uesugi, M., Sekiya, M., & Nakamura, T. 2008, Meteo. & Planet. Sci., 43, 717
Yagishita, M., Choi, S., Watanabe, T., Matsuno. J., & Tsuchiyama, A. 2011, Plasma Conference 2011,22P067-P (abstr.)
Yamamoto, T., & Hasegawa, H. 1977, Progress of Theoretical Phys., 58, 816

120

TABLE 1
Parameters for the estimation of surface energy
surface energy

symbol[J/m^2]

references

iron

γi

1.75

Dyson 1963

chondruleγc0.4King 1951

boundaryγb

1.9

sample 3 (Uesugi+ 2008)

grain radius

	
 

[nm]

	
 

silicaters

50

typical

iron

ri

20

typicalTABLE 2
Λ parameters for astronomical environments and the ITP systemnovae
AGB stars

c1(0) [cm-3]

τT [sec]Λ

5×104 – 5×106

9×106

0.9 – 9

10

11

10
12

presolar nebula (at 0.1 AU)

7×10

this study (ITP)

2.1×10 16

7

3×10

10 6
8

0.12

3×10 94×10 3

TABLE 3
Parameters used for estimating Λ value in the ITP system	 Symbolnon dimensional parameterΛtimescale of saturation

τsat

sec

2.3×103timescale of collision

τcoll

sec

6.1×107timescale of cooling

τT

sec

0.12

equiribulium temperature

Te

K

1000

h: latent heat of silicate

h/kTe-1

sec

50

mass of monomermg

1.7×10-22

radius of monomera0ÅDimension

Value
4×103

2.3
3

initial monomer concentration

c1(0)

/m

mean thrmal velocity

<v>

m/s

sticking probability

αsBoltzman cofficientk

2.1×1022
115
1

m2 kg /s2 /K

1.38×10-23

121

Fig. 1: Schematic illustration of the induction
thermal plasma (ITP). The starting material was
taken into the feeder and carried by Ar gas. Ar(60 L/min)-He (5 L/min) gas was used for sheathgas. 5 L/min of Ar gas was used for inner gas.
Plasma torch was made at 1 bar. Radio frequency(RF; 4 MHz) of 33 kW was supplied. Run
products adhering on the wall of the innerchamber were collected.

Fig. 2: A XRD pattern of the run product. They
were exposed to CuKα radiation (λ=1.54056 Å ).
Two peaks were identified with α-iron and a halo
was caused by amorphous silicate.

Fig. 3: Bright-field micrographs of the run
products by using TEM. (a,b) stereogram of therun product. (c) single grain composed ofamorphous silicate embedded with an iron grain.

122

Fig. 4: Size distributions obtained by using a
TEM image (solid lines; for ~100 grains ≥ 30 nm
of silicate and ≥ 10 nm of iron) and calculated by
Yamamoto and Hasegawa (1977) (dashed lines).
No obvious peaks were observed in the run
product however sharp peaks were shown in thecalculation results.

(IN)

Fig. 5: Bulk and the individual chemical
abundance of the run product plotted in
Si-Mg-Fe ternary diagram obtained by using
the EDS equipped with the FE-SEM and theTEM, respectively (atomic ratio). px: pyroxene,
ol: olivine.

(OUT)

Fig. 6: Schematic drawing of the states of an iron grain and a silicate grain. (IN) a iron grain is insidea silicate grain. (OUT) a iron grain is outside a silicate grain. γs, γi, andγb are the surface energyof silicate, iron, and the boundary of the silicate and the iron, respectively. rs, r’s, and ri are the radiusof silicate for the “IN” state, silicate for the “OUT” state, and the iron, respectively. Et,IN: total surface
energy for the “IN” state, Et,OUT: total surface energy for the “OUT” state.

123

Coulomb-Crystallization and Behavior of Fine Particles inExperimental Plasmas
Yasuaki HAYASHI

Department of Electronics and Information Science, Kyoto Institute of Technology,
Matsugasaki, Sakyo-ku, Kyoto 606-8585, Japan

Because the behavior of fine particles in plasmas on the ground can be directly
observed, the experimental results may more or less bring useful information to
the study on grains in space, although conditions of plasmas are different to eachother.
Mono-disperse spherical fine particles can form Coulomb crystal when they
are included in great abundance in a glow-discharge plasma. Figure 1 shows apicture of a Coulomb crystal formed with spherical carbon fine-particles of 1.4
micron in diameter in methane plasma [1]. Such crystal-like arrangement offine-particles is formed by the non-local coagulation force of positive ionssurrounding negatively charged fine-particles.
Because the theory ofstrongly-coupled plasma, which includes Coulomb crystal formation, depends on
dimensionless parameters, it can be applied to a lot of phenomena related with
fine-particles and plasma in space.

Fig.1 3-dimensional Coulomb crystal formed by fine-particle plasma [1].
The structure of the crystal is body-centered tetragonal.

Behaviors of individual fine particle in glow discharge plasmas have been
investigated by the analyses of video images. We investigated the behavior offine-particles in planar magnetron plasma. Figure 2 shows a top view of
divinylbenzene fine-particles of 2.27 microns in diameter in helium plasma
revolving to the direction of E×B drift [2]. Under a pressure of 100 Pa and a124discharge power of 2 W, fine particles were observed to revolve with the speed of
1.3 cm/s in the direction of E×B, forming particle strings in the vertical direction in
the sheath. A velocity of revolution of a fine particle was calculated under theeffects of drag force by E×B drift ions and friction force by neutral molecules to be
42 cm/s, which is around 30 times larger than the speed observed experimentally.
The Hall parameter of ions is 0.02 under the conditions of the experiment. Thus,
the revolution of fine-particles to the direction of E×B drift should be caused by thedrag force of streams of slightly magnetized ions.

Fig.2 Top view of revolution of 6.5-micron fine particles by E×B drift in helium plasma
under pressure of 100 Pa and RF power of 2 W [2]. Video frames were piled up for 0.1 second

Fine particles must have rotational angular momentum in plasma, which may
more or less influence their behavior. However, the value of rotational angularmomentum has not been evaluated. We evaluated the magnetic moment of arotating fine particle by the analysis of video image. Disc-shaped fine-particles, ofwhich the average diameter was 6 microns, were injected into an RF plasma ofargon under a discharge power of 2 W and a pressure of 100 Pa. Polarized laser
light of 532 nm in wavelength was irradiated. The scattered laser light fromdisc-shaped fine particles at a right angle was observed using a high-speed CCD
video camera. The speed of rotation was determined from the period of intensitychange of scattered light to be 0.03 sec/rotation as shown in Fig.3. When theparticle is charged by -10,000 e and the electrons are distributed around the
periphery of the disc-shaped fine particle, the magnetic moment is calculated to be
3 × 10-25 A・m2.

Fig.3

Change of intensity of light scattered from fine-particles.

[1] Y.Hayashi: Phys. Rev. Lett. 83, 4764 (1999).
[2] Y. Hayashi, Y. Mizobata and K. Takahashi: J. Plasma Res. SERIES, 8, 298 (2009).

125

NUMERICAL SIMULATIONS OF COLLISIONS OF DUST AGGREGATESCOMPOSED OF PARTICLES WITH SIZE DISTRIBUTIONS
Koji Wada 1, Hidekazu Tanaka 2, Toru Suyama 3, Hiroshi Kimura 4, and Tetsuo Yamamoto2,4

Abstract: We carry out numerical simulations of collisions of dust aggregates composedof particles with size distributions to investigate the collisional growth condition of dust in
protoplanetary disks. In this preliminary study we find that the particle size distributions
tend to encourage dust growth, but we need to further study to confirm such a finding.

Introduction:

In protoplanetary disks planetesimals are believed to be formed from dust aggregatesconsisting of submicron grains, but their growth process is unclear so far. One of the main problems inplanetesimal formation is the feasibility of dust growth through collisions at velocities up to several tensof m/s [e.g., 1]. Recently we performed numerical simulations of dust aggregate collisions and found that
dust aggregates can grow through high-velocity collisions even at 50 m/s if they consist of ice particles[2]. However, silicate dust aggregates are still difficult to grow through collisions since their criticalcollision velocity for collisional growth is only several m/s. We need to seek some kind of promisingmechanism to enable silicate dust growth. One of factors we have not yet taken into account in our
numerical simulation is the size distribution of monomer particles. According to a particle-interaction
model, so-called JKR theory [3], the critical collision velocity for sticking of two monomers is
proportional to R-5/6 with R being their reduced radius, suggestive of an increase in the growth possibilityof dust aggregates consisting of smaller particles. In fact, interstellar extinction observation suggests a
power-law size distribution of interstellar dust grains with an exponent of -3.3 to -3.6 in a range of 0.025– 0.25 micron radius [4]. Such a particle size distribution may contribute to enhancement of collisional
growth of dust aggregates. On the other hand, energy to break up a pair of particles in contact is
proportional to R4/3. This means connections between small particles are weaker than those between large
ones. In that sense, a power-law size distribution of constituent particles may have a negative effect forcollisional dust growth. In this study, we carry out numerical simulations of aggregate collisions toinvestigate the influence of particle size distributions on the collisional growth of dust aggregates. Basedon the numerical results, we will discuss the possibility of planetesimal formation through direct
collisions of dust aggregates in protoplanetary disks. Preliminary results indicate that particle size
distributions tend to encourage dust growth.

1Planetary Exploration Research Center, Chiba Institute of Technology (wada@perc.it-chiba.ac.jp)
Institute of Low Temperature Science, Hokkaido Univerisity
3
Nagano City Museum
4Center for Planetary Sciecne2

126

Numerical model and settings: We performed 3D simulations of aggregate collisions by the use of
the particle interaction model and the numerical code developed in the previous papers [2,5,6]. We
directly calculate the motion of each particle, taking into account all mechanical interactions between
particles in contact. The contact theory of adhesive elastic spheres determines the interactions for each
degree of motion (normal motion, sliding, rolling, and twisting). Energy dissipates at the moments ofcontact and separation of particles because of the excitation of elastic waves. When the displacements due
to sliding, rolling, and twisting exceed their elastic limits, the mechanical energy is also dissipated. The
amount of energy dissipation is proportional to the critical displacements.
In this study, as a first step, we assume a power-law size distribution of constituent particles:𝑛(𝑟)𝑑𝑟 ∝ 𝑟 −3.5 𝑑𝑟,

where n(r)dr is the number of particles in the range of radius r to r+dr. The power-law index of -3.5 is inaccordance with the values obtained from interstellar dust observation [4]. Along with this size
distribution, we prepare aggregates composed of particles with a bimodal distribution and a continuous
distribution. For a bimodal distribution, we mix large particles of 0.1µm radius and small ones of 0.025µm radius with a number ratio 1/32. Then we prepare two types of aggregates composed of 3300 particles(100 large particles and 3200 small ones) and 6600 particles (200 large and 6400 small) by means ofrandom sticking of the particles one by one, namely ballistic particle-cluster aggregation (BPCA) process.
Their masses are equivalent to 150 and 300 mu, respectively, where mu is the mass unit defined by themass of 0.1µm-radius particle. On the other hand, for a continuous distribution, we set the largest radiusof 0.2 µm and the smallest of 0.025 µm and determine particle radii along with the power-law distribution.
Then, we prepare an aggregate composed of 6600 particles by means of BPCA process. This aggregate
has an equivalent mass of 201.6 mu. As a result, we prepare three kinds of aggregates, two of which are
shown in Figure 1. They are collided with the same type aggregates in the head-on direction and theirnumerical results are compared with our previous ones [2] obtained from collisions of aggregatescomposed of equal-sized particles of 0.1 µm radius. All particles in this study are assumed to be ice

Figure 1: Examples of initial aggregates prepared in this study. Left panel shows aggregates composed of 200
particles of 0.1µm radius and 6400 particles of 0.025µm radius (bimodal distribution), while right panel shows
aggregates composed of 6600 particles with a continuous power-law size distribution (the range of particle radius
is 0.025 – 0.2 µm, and the power index is -3.5).

127(Young’s modulus 7 GPa; Poisson’s ratio 0.25; material density 1000 kg/m3; surface energy 100 mJ/m2).
Since we examine the growth possibility of dust aggregates at high-velocity collisions, we simulate
aggregate collisions at collision velocities ucoll in the range of 20 – 700 m/s. We can estimate the results
for silicate aggregates from those for ice by replacing ucoll by one order of magnitude lower value [2].

Results and discussion: Figure 2 shows some examples of snapshots of moderate- and high-velocity
collisions of our prepared aggregates. Large parts of resultant aggregates are likely to survive even at highcollision velocities ucoll ~ 50 and 200 m/s. Fragments produced at these collisions are composed ofrelatively large particles, a particular feature of aggregate collisions with particle size distributions.
To examine the collisional growth possibility of aggregates, we introduce the growth efficiency fdefined by
𝑓 = �𝑀large − 𝑀init �⁄𝑀init ,where Mlarge is the mass of the largest fragment at a collision and Minit is the mass of an initial aggregate.
The growth efficiency counts how much mass, normalized by Minit, are gained (or lost, if the value isnegative) for an initial aggregate through a collision. Figure 3 shows f plotted as a function of ucoll for thethree types of aggregates, together with the previous results [2] obtained from head-on collisions of
BPCA clusters that are composed of ice particles with the same radius of 0.1µm and have equivalentmasses of 500, 2000, or 8000 mu. First, compared with the previous BPCA clusters, the growth efficiency

Figure 2: Snapshots of collided aggregates shown in Figure 1. Upper two panels are of the bimodal distribution
at ucoll =55.3 m/s (left) and 180 m/s (right) and lower two are of the continuous distribution at ucoll = 49.3 m/s(left) and 197 m/s (right).

128

f for the aggregates with particle size distributions in this study becomes large for the same ucoll eventhough their small masses. This indicates that aggregates with particle size distributions tend to be hardfor disruption and that small particles result in positively working for collisional growth of dust
aggregates. On the other hand, we do not see significant difference in f between aggregates with thebimodal particle size distribution and the continuous one. In addition, f seems to be independent of theaggregate mass for the cases with particle size distributions, while aggregates composed of equal-sizedparticles have a mass-dependence in f. As a result, the critical collision velocity ucoll,crit, which is given by
the collision velocity when f = 0, is found to be ~ 200 m/s for the ice aggregates with particle size
distributions used in this study. Taking into account the difference in ucoll,crit between ice and silicate
aggregates, ucoll,crit for silicate aggregates is estimated to be ~ 20 m/s. This critical collision velocity would
be high enough to survive in optically thick and cold protoplanetary disks [7], but it should be cautionedthat this result is obtained from head-on collisions: ucoll,crit will be reduced if offset collisions are taken
into account. We need to further study in a wide range of the aggregate mass, the particle size distribution,and the collision direction (impact parameter) to investigate whether or not particle size distributions
enhance the collisional growth of dust aggregates.
Furthermore, the size distribution of fragments at collisional disruption of dust aggregates is expected

Figure 3: Growth efficiency f vs. collision velocity obtained from this study (for aggregates composed ofparticles with size distributions, represented by solid lines) and from the previous study (for aggregates composedof equal-sized particles, represented by dashed lines). Numbers in each parenthesis denote masses of colliding
aggregates in mass unit (mu): for example, “150+150mu” means two aggregates of 150 mu collide with eachother.

129

to be affected by the size distribution of constituent particles as can be seen in Figure 2. Since fragment
size distribution is an important factor in observing protoplanetary or debris disks, we will investigate thefragment size distributions obtained from our simulations in the future.

References:[1] Weidenschilling, S. J., & Cuzzi, J. N. (1993), in Protostars and Planets III, 1031.
[2] Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. (2009), ApJ 702, 1490.
[3] Johnson, K. L. (1987), Contact Mechanics.
[4] Mathis, J. S., Rumpl, W., & Nordsieck, K. H. (1977), ApJ 217, 425.
[5] Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. (2007), ApJ 661, 320.
[6] Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. (2008), ApJ 677, 1296.
[7] Kusaka, T., Nakano, T., & Hayashi, C. (1970), Prog. Theor. Phys. 44, 1580.

130

Ionization Degrees with Dust Grains in
Circumplanetary Disks
Yuri I. Fujii1 , Satoshi Okuzumi1 , and Shu-ichiro
Inutsuka1
1

Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya,
Aichi 464-8602, Japan
E-mail: yuri.f@nagoya-u.jp
Abstract. During the formation of gas giants, disks appear around them.
They are called circumplanetary disks. The most promising mechanism ofaccretion is the turbulence driven by magnetorotational instability (MRI).
MRI occurs when the ionization degree is high enough for magnetic fields
to couple to gas. We calculate the ionization degrees to investigate if MRI
really causes accretion of circumplanetary disks. Our result suggests that theionization degree in circumplanetary disks is too low for MRI to be active.

1. Introduction
There are only a few studies on the ionization degrees in circumplanetary disks,
although a number of papers concerning ionization degrees in protoplanetarydisks have appeared in the literature [1] [2] [3]. However, it is important tounderstand the structure and evolution of circumplanetary disks, because they
are thought to be the sites of satellites formation.
The magnetohydrodynamic (MHD) turbulence driven by the magnetorotational instability (MRI) is, at the present, thought to be the most promisingmechanism for angular momentum transport in protoplanetary disks. However,
we are not sure if the MRI also plays a central role in accretion in circumplanetary disks. In order to understand disk evolution, we must clarify which part ofthe disk is magnetorotationally unstable.
In this work, we develop a fast and accurate time-dependent calculation
method for the ionization degree in protoplanetary and circumplanetary disks.
Using our calculations of ionization degree, we estimate the size of the
MRI-inactive region. We ﬁnd that almost all regions in popular model ofcircumplanetary disks are MRI-inactive.
2. Developing the calculation method
2.1. Reaction Equations
Gas in protoplanetary and circumplanetary disks is mostly in the form ofneutral hydrogen molecules. However, they are thought to be weakly ionizedby ionization sources such as cosmic rays, X-rays, and ultraviolet radiation. The
resultant ionized particles make secondary ions and molecules, which lead tofurther complex reactions.

2

Fujii, Okuzumi, & Inutsuka

The ionization degree is determined by the balance of chemical reactions,such as ionization, charge transfer, and dissociative and radiative recombination.
Since molecular ions transfer their charge to metal ions rapidly, we need not pay
much attention to reactions of molecular ions in the dense gas. We describe
major heavy metal ions (e.g. Mg+ , Fe+ ) as M+ for simplicity [4].
The eﬀect of dust grains should also be considered. Ions and free electrons arecaptured by dust grains, lowering the ionization degree. With suﬃcient amountof dust, the ionization degree is determined mainly by the electron capture rateof dust grains. We show simpliﬁed reaction equations as follows:dnd (Z)
dt

dnM+
= ζnn
dt

M+ nM+ ne

dne
= ζnn
dt

M+ nM+ ne=∑

kM+ d (Z)nd (Z)nM+ ,

(1)

Z

∑

ked (Z)nd (Z)ne ,(2)

Z

kM+ d (Z)nd (Z)nM++kM+ d (Zked (Z)nd (Z)ne1)nM+ + ked (Z + 1)nd (Z + 1)ne ,1)nd (Z(3)

where nj is the number density of each particle (j = M+ : heavy metal ions, e :electrons, d : dust grains), Z is the charge of dust grains, ζ is the ionization
rate, and M+ and kjd (Z) are the rate coeﬃcients for radiative recombinationand capture of gaseous particles by dust grains, respectively.
Numerical calculation of Equations (1) (3) is not an easy task. As the
maximum value of |Z| becomes greater, the number of terms and equations
increases, and the system of equations becomes more complicated.
2.2. Charge distribution model of dust grains
Since the charge distribution of dust grains can be approximated by a Gaussian
distribution [1], we adopt the following equation for use in Equations (1) (3).
Ndnd (Z) = √exp
2hδZ 2 i[]

(Z hZi)2,2hδZ 2 i(4)where h i shows the averaged value weighed by nd (Z), and Nd is the total number
density of dust grains. We can get the following equations by taking the averageof the last term of Equations (1) and (2), and the ﬁrst and second moments ofEquation (3):
dnM+
dt
dne
dt
dhZi
dt
dhδZ 2 idt= ζnn

M+ nM+ ne

hkM+ d iNd nM+ ,

(5)

= ζnn

M+ nM+ ne

hked iNd ne ,(6)

= hkM+ d inM+

hked ine ,(7)= (hkM+ d i + 2hkM+ d δZi) nM++ (hked i

2hked δZi) ne .

(8)

131

Ionization Degrees with Dust Grains in Circumplanetary Disks

3

The details of our calculations are described in our paper [5]. We have conﬁrmedthat the results of our equations agree with those of the basic equations very well.
3. Ionization Degrees in Circumplanetary Disks
In the context of satellite formation, understanding the evolution of
circumplanetary disks is important. Since the ionization degree is one of
the important factors in disk evolution, we calculate the ionization degree incircumplanetary disks including the charged particle capture by dust grains.
We use an actively supplied gaseous accretion disk model [6] [7]. The
temperature is(

Tcir ' 160

Mp
MJ)1/2 (

r
20RJ

)−3/4

K,

(9)

where Mp is the mass of the central planet, MJ is that of Jupiter, r is the orbital
radius from the central planet, and RJ is the radius of Jupiter. For this model,
we adopt 5  106 yr as the accretion timescale for gaseous disks as heavy as thecentral planet, and assume the disk is vertically isothermal. The surface density
is given as follows with viscosity coeﬃcient = 5  10−3 :(MpΣcir ' 100
MJ

)(

r
20RJ

)−3/4

g cm−2 .

(10)

We assume Mp = MJ . Since circumplanetary disks are located at the orbital
radius of gas giant planet which is distant from the star, stellar X-rays aregeometrically attenuated there. Furthermore, the scale height of circumplanetarydisks is far smaller than that of protoplanetary disks. So, it is diﬃcult for
the X-rays to directly reach circumplanetary disks. When we think aboutcircumplanetary disks, we consider only cosmic ray ionization. The ionization
rate is written as follows:
ζ '

ζCR
2

{

[

exp

]

χ(r, z)
+ exp
χCR

[

Σ(r)

χ(r, z)
χCR

]}

,

(11)

where ζCR  1.0  10−17 s−1 is the ionization rate by cosmic rays in interstellar
space, and χCR  96 g cm−2 is the attenuation length of the ionization rate bycosmic rays. The vertical column density of the gas measured from the outsideof the disk is the following [8]:χ(r, z) =∫ ∞z(r, z)dz [g cm−2 ].

(12)

In this study, we investigate the extent of the MRI-inactive regions in thecircumplanetary disks. We parameterize the vertical component of plasma
beta (at the mid-plane), the ratio of gas pressure and magnetic pressure,β = Pgas /Pmag , the radius of dust grains a, and the dust-to-gas mass ratio
fdg = d /g . We deﬁne the magnetic Reynolds number Rem as the following:
Rem =

2
vAz
,
ηΩK

(13)

132

4

Fujii, Okuzumi, & Inutsuka
4

z/H

3
2
10−3
0.1
1

1
0
5

10

15

20

25

30

r [RJ]

Figure 1. Boundaries of the MRI-inactive region in circumplanetary disks forthe case β = 104 and fdg = 10−2 . The horizontal axis denotes the orbital radiusfrom the central planet, and the vertical axis is vertical extent of the disk that isnormalized by the scale height of the corresponding radius. Contour lines show
the boundaries of the MRI-active region and the MRI-inactive region for eachmodel of dust grain radius a The circumplanetary disks are not expected to bemagnetically active.
where vAz is the vertical component of the Alfvén velocity. If Rem < 1, MRI
does not happen. Figure 1 shows an example of the extent of the MRI-inactive
regions.
For a wide range of model parameters, we ﬁnd that Rem falls below unity
in the entire parts of the circumplanetary disks. Our results suggest that theionization degree of circumplanetary disks is too low for MRI to be active.
References[1] Okuzumi, S. 2009, Astrophys. J., 698, 1122
[2] Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T., 2000, Astrophys. J., 543, 486[3] Ilgner, M., & Nelson, R. P., 2006, A & A, 445, 205[4] Oppenheimer, M., & Dalgarno, A., 1974, Astrophys. J., 192, 29[5] Fujii, Y. I., Okuzumi, S. & Inutsuka, S., 2011, Astrophys. J., 743, 53[6] Canup, R. M., & Ward, W. R., 2002, Astron. J., 124, 3404
[7] Sasaki, T., Stewart, G. R., & Ida, S., 2010 Astrophys. J., 714, 1052
[8] Umebayashi, T., & Nakano, T., 1981, PASJ, 33, 617

133

134

Hydrogen ordering in ice observed from neutron diffraction: Application toplanetary science
Masashi Arakawa
Department of Chemistry, Faculty of Sciences, Kyushu University

Laboratory experiments and infrared observations suggest that a lot of water ice
exists in a crystalline phase in our solar system. At ambient pressure, crystalline ice has
two kinds of structure, ice Ih and ice XI. Ice Ih is normal ice with disordered hydrogenatoms. In contrast, ice XI is a hydrogen-ordered phase of ice Ih, and the orderedarrangement of hydrogen atoms makes ice XI ferroelectric. Existence ofhydrogen-ordered ice in space is the subject of continuing astronomical debate.
Electrostatic forces, caused by the ferroelectricity, increase the sticking probability oficy grains, and might play an important role in grain evolution and planetary formationin space.
To discuss the existence of hydrogen-ordered ice in space, we need toinvestigate the kinetics and formation process of ice XI. For precise understanding of
nucleation and growth process of ice XI, we performed neutron diffraction
measurements of ice doped with various kinds of dopant as a catalyst for thetransformation to ice XI. The neutron diffraction measurements were performed at thehigh-resolution powder diffractometer (HRPD), and the wide-angle neutron
diffractometer (WAND). The HRPD is installed in the research reactor, JRR-3, JAEA,
Japan. The WAND is installed in the research reactor, High Flux Isotope Reactor,
ORNL, USA. All the neutron diffraction measurements were performed on 99.95%
deuterium-substituted ice. Rietveld analysis was carried out for all the diffraction
profiles to obtain arrangement of hydrogen (structure parameter). We used a two-phasemodel, which include ice Ih and XI, and obtained mass fraction of ice XI (f).
Neutron powder diffraction of 0.013 M KOD-doped D2O ice at HRPD. The
doped ice Ih transformed to ice XI after annealing at 57 and subsequently at 68 K. The f
value of the doped ice, which had once experienced being ice XI (f = 0.23), was larger
than that of the doped ice, which had never experienced being ice XI (f = 0.14). Results
indicate that small hydrogen-ordered domains remained in the ice Ih, which had oncetransformed to ice XI, and accelerated the phase transition from ice Ih to ice XI.

135

To understand this phenomenon in detail, we performed time-resolved neutron
diffraction measurements on 0.1 M NaOD-doped ices, and investigated whether
formation of ice XI occurs at temperatures above 65 K, which is out of the nucleation
temperatures. Formation of ice XI occurred at 70 and 72 K when the ice had experience
being ice XI in the past. This result suggests the presence of a template acting as thenuclei of ice XI. Thus, small domains with ordered hydrogen, which cannot be detectedwith neutron diffraction, can exist above the boundary temperature between ice Ih and
XI. The small ordered-domain is called “memory” of hydrogen ordered structure
because of the residual structure of ice XI. In the previous studies, ferroelectrichydrogen-ordered ice is considered to exist in the very narrow temperature range, but
our experiments show a possibility that the hydrogen ordered ice exists in wider area inspace.
In this presentation, we also discussed hydrogen ordering in ice at high pressureusing neutron diffraction experiments. Our results indicate that hydrogen-orderedstructure of high-pressure ice VI is ferroelectric. Hydrogen-ordered ferroelectric icesmight exist not only in icy grains and icy bodies’ surface, but also in icy bodies’interior.

136

Proceeding of the 29th Grain Formation Workshop in CPS center (2011/11/9-11)

Evaporation of Icy Planetesimals due to Planetesimal Bow Shocks
in a Protoplanetary Disk
Kyoko K. Tanaka1 , Tetsuo Yamamoto1 , Hitoshi Miura2 , Makiko Nagasawa3 ,Taishi Nakamoto3 , and Hidekazu Tanaka1
1 Institute of Low Temperature Science, Hokkaido University, Sapporo 060-0819, Japan
2 Department of Earth Sciences, Tohoku University, Sendai 980-8578, Japan
3 Earth and Planetary Sciences, Tokyo Institute of Technology, Tokyo 152-8551, Japan

In a protoplanetary disk, planetesimals grow to planets by mutual collisions and accumulations.
The gravitational interactions among the planetesimals increase eccentricities of the planetesimal
orbits. When a relative velocity between the disk gas and planetesimal exceeds a sound velocity of
the gas, a bow shock wave is produced on the leading side of the planetesimal.
It has been shown that heating by the planetesimal bow shocks in the nebular gas played a keyrole in formation of particles found in meteorites and interplanetary dust. For instance, the shockheating leads to formation of chondrules, millimeter-scale igneous silicate spheres in chondrites,by melting dust in the protoplanetary disk (Hood 1998; Weidenschlling et al. 1998; Ciesla et al.
2004; Hood et al. 2009). Furthermore, it makes a part of ﬁne dust particles evaporate (Miuraand Nakamoto 2002, 2005, 2006) and re-condesation of evaporated vapor forms a lot of small dust.
Miura et al. (2010) showed that it is possible by a planetesimal bow shock to form various types of
cosmic crystals, which are ﬁne silicate crystals observed in chondritic meteorites and interplanetary
dust particles (IDPs).
So far, the attention has mainly been paid on thermal evolutions of dust particles in theshocked region in the previous studies. However, we note a possibility that the planetesimal bow
shock leads to heating and evaporation of the planetesimal itself. A similar process is found in
ablation of the planetesimals penetrating through the atmosphere of a protoplanet (e.g., Podolak
et al. 1988; Pollack et al. 1986). Heating and resultant evaporation by the planetesimal bow shocks
will suppress growth of planetesimals. Furthermore, cooling of the vapors thus produced will form
small dust particles by re-condensation. Those dust clouds in the disk may be found in the infrared
spectra of the protoplanetary disks.
In the present study, we focus on the planetesimal heating and evaporation by the planetesimal
bow shocks. Figure 1 shows an illustration of the planetesimal evaporation by the bow shock wave.
When the planetesimal velocity relative to the disk gas becomes supersonic, a bow shock forms infront of the planetesimal. Due to the shock wave, the gas is heated at the shock front. The gas

137

–2–
blows onto the planetesimal surface. The surface material of the planetesimal heats due to theheat transfer from the gas. The surface material evaporates when the temperature rises over theevaporation temperature. The evaporated vapor re-condenses and forms small dust particles. As
a result of the evaporation, the planetesimal is shrinking.
We evaluated the surface temperature and evaporation rate of the planetesimal using a simple
model of planetesimal evaporation by the planetesimal bow shock. The time evolution of theplanetesimal mass m is given by
dm
= J,
dt

(1)

where J is the evaporation rate of a planetesimal. In our model, the rate of heat transfer from thegas to the planetesimal surface is given by
Fe = αρ0 vp cp (Ts − Tp ),(2)

where Ts is the gaseous temperature of the post shock region, Tp is the surface temperature of theplanetesimal at the stagnation, ρ0 and cp are the density and speciﬁc heat at constant pressure ofambient gas, and vp is the relative velocity between the ambient gas and the planetesimal. Also α
is a non-dimensional parameter, referred to as the Stanton number, expressing the eﬃciency of thethermal conductivity (van Driest 1958; Schlichting 1979).
We show one typical example of the calculations in Figure 2 which shows the surface temperature of the planetesimal at the stagnation point Tp , the gaseous temperature of the post shock
region Ts , and the ambient radiation temperature Tdisk (= 75K). In Figure 2, we set α = 0.1 and
0.01. As shown in Figure 2, the temperature of the planetesimal begins to increase when the relative
velocity vp exceeds 0.6kms−1 which is the sound velocity of the gas. Although the temperatures ofboth Ts and Tp increase with increasing of vp , the elevation of Tp is suppressed in the vicinity of−1 because of the cooling by latent heat release.
200K for vp >∼ 5kms
We applied the model of the planetesimal evaporation to the formation stage of protoplanets.
According to the theory of formation of planets, the terrestrial planets and cores of the gas giant
planets are formed through the accretion of planetesimals (e.g. Hayashi et al. 1985). According to
the standard scenario (Wetherill and Stewart 1989), growth of a planet pass through a so-calledrunaway growth stage. In this stage, the large planetesimals grow faster than small ones; therate of mass growth of each object is an increasing function of its own mass. As a result, theplanetary embryos grow to protoplanets accumulating ambient planetesimals with keeping theirorbital separations (Kokubo and Ida 1998, 2002). In this oligarchic growth stage, the eccentricities
and inclinations of the surrounding planetesimals are pumped up by the gravitational interaction
with the protoplanet. On the other hand, the eccentricities and inclinations of the planetesimal
dumps due to the gas drag (Tanaka and Ida 1999). Balancing the stirring and damping rates, theequilibrium value of e is obtained (Tanaka and Ida 1999; Kobayashi and Tanaka 2010).
Figure 3 shows the evaporation time of the planetesimal of 100km size as the function of the

138

–3–
semi-major axis a, where we deﬁne the evaporation time of the planetesimal tevap as(tevap ≡

1 dm
m dt

)−1
=

m
.
J

(3)

In Figure 3, the mass of the protoplanet is assumed to be ten Earth mass. The evaporation times are
below 107 yr from the snow line (0.5AU) to 4 AU for α = 0.1. The evaporation time is shorter withthe decrease of the semi-major axis. The disk region where the evaporation of the planetesimals is
eﬀective becomes narrower for smaller value of α.
Our results show that the icy planetesimals evaporate eﬃciently in the planetary oligarchic
stage of the protoplanets, where strong shocks are generated by the gravitational perturbations fromthe protoplanets. There are several important implications to planetary formation and observation.
The growth of a protoplanet is suppressed owing to the insuﬃcient accretion of planetesimals onto
the protoplanet. The result suggests that the evaporation of the planetesimals disturbs the growthof the Jupiter type planet core. It was shown to disturb the growth of the Jupiter type planet coreby the eﬀect of destruction because the collision among the planetesimals leads the destruction of
the planetesimals (Kobayashi et al. 2010). Our results suggest that the evaporation as well as thedestruction of the planetesimal would be eﬀective, especially for the inner disk region.
The evaporated vapor re-condenses and forms small dusts. Moreover, solid materials would
release with evaporating gas because dust particles are scattered by the strong gaseous ﬂow. This
process gives an inﬂuence on the observations of the protoplanetary disks. The observations of the
protoplanetary disks by SED (spectral energy distribution) have shown to remain dust rich for up
to several million years, although the theoretical modelings show that dusts growth proceeds invery short time (Birnstiel et al. 2009). The dust clouds formed by the planetesimal bow shocks
may explain the observational results of the SED.
In the present study, we focused on the growth stage of the protoplanets where the surrounding
planetesimals around a protoplanet have high relative velocities. The stage after the formation ofgas giant planets is an other stage increasing eccentricities and inclinations of the planetesimals.
The secular resonances due to the gas giant planets give strong perturbations on the asteroid region(Nagasawa et al. 2000). In such a stage, the evaporation of the planetesimal would occur more
eﬀectively.

1.

REFERENCES

Birnstiel, T., Dullemond, C. P., & Brauer, F., 2009, A&A 53, L5
Ciesla, F. J., Hood, L. L., & Weidenschilling, S. J. 2004, Meteorit. Planet. Sci., 39, 1809
Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, Formation of the solar system. In Protostars
and Planets II, eds. D.C. Black and M.S. Matthews, Univ. of Arizona Press, 1100-1153.
Hood, L. L. 1998, Meteorit. Planet. Sci., 33, 97

139

–4–
Hood, L. L., Ciesla, F. J., Artemieva, N. A., Marzari, F., & Weidenschilling, S. J. 2009, Meteorit.
Planet. Sci., 44, 327
Kobayashi,H., & Tanaka, 2010, Icarus, 206, 735
Kobayashi,H., Tanaka, H., Krivov, A. V., & Inaba S., 2010, Icarus, 209, 836
Kokubo, E. & Ida, S. 1998, Icarus, 131, 171
Kokubo, E. & Ida,S. 2002, ApJ, 581, 666
Miura H., Nakamoto T, and H. Susa, 2002, Icarus 160, 258
Miura H., and Nakamoto T. 2005, Icarus 175, 289
Miura H., and Nakamoto T. 2006, ApJ, 651, 1272
Miura, H., K. K. Tanaka, T. Yamamoto, T. Nakamoto, J. Yamada, K. Tsukamoto, and J. Nozawa,
ApJ., 2010, 719, 642
Nagasawa, M., Tanaka, H., & Ida, S. ApJ, 2000, 119, 1480
Pollack, J. B., Podolak, M., Bodenheimer, P., & Christoﬀerson, B. 1986, Icarus, 67, 409
Podolak, M., Pollack, J., & Reynolds, R. T., 1988, Icarus, 73, 163
Schlichting, H. 1979, in Boundary-Layer Theory (seventh edition), McGraw-Hill, Inc. 5
Tanaka H., and Ida S. 1999, Icarus,139, 350
Weidenschilling, S. J., Marzari, F., Hood, L. L. 1998, Science, 279, 681
Wetherill, G. W., & Stewart G. R. 1989. Icarus, 77, 330
van Driest, E. R. 1958, Zeitschrift fr Angewandte Mathematik und Physik (ZAMP), 9, 233

140

–5–

evaporation fromplanetesimal surface
re-condensation
disk gas ﬂow

icy
ic
cy
p
planetesimallaanetesimmalboshrinkage

w

sh

oc

k

Fig. 1.— A illustration of a planetesimal evaporation due to planetesimal bow shock. A shock wave
is produced on the leading side of the planetesimal when a relative velocity between the disk gas andplanetesimal exceeds a sound velocity of the gas. The planetesimal surface heats and evaporates.
The evaporated vapor re-condenses and forms small dusts particles.

141

–6–

3

[K]

10

α=0.1
10

α=0.01

2

2

4

6
-1

[km s ]
Fig. 2.— The temperature of the planetesimal surface at the stagnation point Tp and gaseous
temperature of the post shock region Ts , where we set the Stanton number α = 0.1 (solid curve)and 0.01 (dotted curve). The disk temperature Tdisk is 75 K and gaseous density is ρ0 = 1.0 ×
10−10 gcm−3 .

142

10

8

10

6

10

4

0

α=0.01

Snow line

tevap [yr]

–7–

α=0.1

mp=10 ME
R=100kmAsteroid region

2

4a [AU]

Fig. 3.— The evaporation time of the planetesimal of 100 km radius as a function of the semi-major
axis a. The mass of the protoplanet is set to be ten Earth mass.

143

Structure of Wind-Wind Collision Shocks in Massive Binaries
Atsuo T. Okazaki
Faculty of Engineering, Hokkai-Gakuen University, Toyohira-ku, Sapporo 062-8605, Japan

Abstract
In binaries consisting of two massive stars, strong stellar winds from these stars collide
each other. Shocked regions produced by the wind-wind collision is thought to be a place
where not only the particle acceleration and resulting high-energy emission but also the dust
formation takes place. Studying the wind-wind collision shocks is, therefore, of multi-facet
importance. Numerical simulations suit it very well, given that massive binaries, in general,
have large orbital eccentricities. In this paper, we report on the structure of the wind-windcollision shocks in two massive binaries, η Carinae and WR 140, obtained from 3-D Smoothed
Particle Hydrodynamics (SPH) simulations.

1.

Introduction

Massive stars have powerful winds driven by strong line radiation. In binaries consisting of two massive
stars, two winds collide and form complicated shock structures, from which high energy emission arise.
Colliding winds occur in binaries with OB stars, LBV (luminous blue variable) stars, and Wolf-Rayetstars. Wind-wind collision also occurs in gamma-ray binaries comprised of a massive star and a pulsar
with relativistic wind..
From 1990s through early 2000s, large progress has been made in dynamical modeling of colliding
winds. Most of those works were two-dimensional, but there were also studies using three dimensionalmodels. With the help of rapid progress of computer power, recent eﬀorts in dynamical modeling ofcolliding winds are mainly made in 3-D [Grid-based simulations: Pittard & Stevens (1999); Walder &Folini (2003); Lemaster et al. (2007); Pittard (2009); Parkin et al. (2011); SPH simulations: Okazaki etal. (2008)]. Although 3-D simulations are sitll expensive, having a 3-D dynamical model is essential tounderstanding of detailed structure of colliding winds in massive binaries and in gamma-ray binaries aswell. Given that several Wolf-Rayet binaries show dust emission originated from wind region, dynamical
modeling of wind-wind collision will also play an important role to understand the origin of dust emission
in these systems.
2.

Basic Geometry of Wind-Wind Collision Shocks

A key parameter that determines the bow shock structure in colliding winds is the ratio of wind
momentum ﬂuxes, η, which is given byη=Ṁ2 v2,Ṁ1 V1(1)where Ṁ1 and v1 are the mass loss rate and velocity of the primary wind, respectively, and Ṁ2 and v2
are those of the secondary wind.
In a simple 2D model that ignores the orbital motion, the ram pressure balance then implies that, fora binary separation D, the interface should be located at a distance
d=

D
√
1+ η

(2)

from the primary (Stevens et al., 1992; Canto et al., 1996). In adiabatic shocks, the interaction surface
approaches a cone with the half opening angle,
(
)1/4
(3)
θ = 2 tan 1 ηmeasured from the primary (Gayley, 2009). The shock is wrapped around the star with the weaker
wind.

144

2

A. T. Okazaki
Table 1. Stellar, wind, and orbital parameters of η Carinae and WR 140

Parameters
Mass (M )Radius (R )Mass loss rate (M yr 1 )Wind velocity (km s 1 )Wind temperature (K)Orbital period Porb (d)Orbital eccentricity eη Carinaeη Car Aη Car B9030
9030
2.5  10 4
10 5
500
3,000
35,00035,0002,024
0.9

WR 140
O4-5V
WC7
50
19
12
13
1.2  10 6 3.8  10 5
3,200
2,860
42,00042,0002,899
0.881

If the thickness of shocked region is ignored, the interaction front with orbital motion becomes an
Archimedean spiral. Such a spiral geometry has been found in several Wolf-Rayet binaries as dust
emission (e.g., Tuthill et al., 2008). They are binaries consisting of a Wolf-Rayet star and an O star.
Thus, colliding winds in a massive binary (WR+O) can form dust, which suggests strong cooling in postshock regions.
In eccentric binaries, the spiral geometry of the interaction region becomes asymmetric, because of thephase-dependent orbital velocity. The asymmetry is more remarkable for a higher orbital eccentricity(Parkin & Pittard, 2008).
3.

3-D Numerical Modeling of Wind-Wind Collision in Highly Eccentric Binaries

3.1.

Numerical model

The simulations presented below were performed with a 3-D SPH code, which is a particle methodthat divides the ﬂuid into a set of discrete “ﬂuid elements” (=particles), and is ﬂexible in setting various
initial conﬁgurations. The code is the same as that used by Okazaki et al. (2008) (see also Bate et al.,
1995; Okazaki et al., 2002). except that the current version takes into account radiative cooling withthe cooling function generated by CLOUDY 90.01 for an optically thin plasma with solar abundances(Ferland, 1996). Moreover, it emulates the empirical β-velocity wind given by
vi = v∞,i (1βRi /r) i(i = 1, 2),(4)where βi is a constant, v∞,i is the terminal velocity of the wind, and Ri is the stellar radius. Here, i = 1,2
denotes the primary and secondary, respectively. Using a variable smoothing length, the SPH equations
with a standard cubic-spline kernel are integrated with an individual time step for each particle. We
adopt standard values of artiﬁcial viscosity parameters, SPH = 1 and βSPH = 2.
3.2.

Supermassive binary η Carinae

η Carinae is one of the most luminous and massive stars in the Galaxy. It has exhibited a series ofmass ejection episodes, the most notable of which was the Great eruption in the 1840’s when the starejected mass of > 10 M , from which the Homunculus nebula formed. Its current mass and luminosity
are M  102 M and L  5  106 L , respectively.
The spectrum of η Car is rich in emission lines and has no photospheric lines. Damineli (1996) founda 5.5 yr periodicity in the variability of the He I 10830Å line. Later, other optical lines were also found toshow variations with the same periodicity. In the X-ray band, the ﬂux exhibits particularly interesting,periodic variability: After a gradual increase toward periastron, the X-ray ﬂux suddenly drops to aminimum, which lasts for two to three months. It then recovers to a level slightly higher than that at
apastron (see Fig. 1 of Corcoran et al., 2010). All these variations are consistent with the system beinga long-period (Porb = 2, 024 days), highly eccentric (e  0.9) binary. The X-ray emission is considered
to arise from the wind-wind collision shocks.
There is also a claim that η Car forms dust in post shock regions. Recently, Smith (2010) analyzed
the IR emission and found that η Car shows broad peaks in JHKL photometry, roughly correlated withtimes of periastron passage and that these peaks have IR SEDs consistent with emission from hot dustat 1400―1700 K. According to Smith (2010), the excess SEDs are inconsistent with the excess being
entirely due to free-free wind or photospheric emission.
Table 1 summarizes the stellar, wind, and orbital parameters adopted in our simulations. With theseparameters, the ratio η of the momentum ﬂuxes of the winds from η Car A and B is η  4.2. The

145

Structure of Wind-Wind Collision Shocks in Massive Binaries

(a)

3

(b)

Fig. 1. Interface geometry of colliding winds of η Car at t ∼ +200 d (200 days after the periastron) within r = 100a:(a) 2-D density maps in the orbital plane (left) and the plane perpendicular to the orbital plane and through the
major axis of the orbit (right) and (b) the 3-D plot of the logarithmic density.

Fig. 2. Density (upper panels) and temperature (lower panels) distribution of η Car’s wind-wind collision region
in the orbital plane in the vicinity of the stars (r ≤ a). The orbital phases are, from left to right, t = −32 d,
t = +1 d, and t = +30 d. Optically-thin radiative cooling is taken into account, with the floor temperature at
35,000 K. Unperturbed wind has β = 1 velocity distribution, which is the value consistent with observed windvelocity distribution for massive stars.

parameters adopted here are consistent with observations, except for the wind temperature of η Car A,
where for simplicity we take the same temperature as that of η Car B wind. Note that the eﬀect ofwind temperature on the dynamics of high-velocity wind collision is negligible.
In Fig. 1, we present a large-scale density structure of the colliding winds. The ﬁgure shows thelogarithmic density distribution of the colliding winds within r = 100a [or 0.7 arcsec ( 3, 000 AU) in
one dimension] of η Carinae. The orbital phase of these plots is 0.1 (about 200 days after the periastronpassage). The two left panels are the 2-D density maps in the orbital plane and the plane perpendicular

146

4

A. T. Okazaki

Fig. 3. Interface geometry of colliding winds of WR 140 at t = +200 d within r = 100a: 2-D density maps in theorbital plane (left) and in the plane perpendicular to the orbital plane and through the major axis of the orbit(right).

to it, while the right panel shows the density distribution in 3-D. From Fig. 1, we note that the lower
density, fast wind from the secondary carves out a large-scale, spiral cavity in the higher density, slow
wind from the primary. Because of high orbital eccentricity, the shape of the cavity is very asymmetric:
It is a thin layer on the periastron side with respect to the primary, whereas it occupies a large volume
on the apastron side.
In the simulation shown in Fig. 1, no radiative cooling is taken into account, so that the post shock
temperature is very high (T  108 K) and no dust formation is possible. On the other hand, Fig. 2
shows snapshots from a simulation where the radiative cooling is taken into account. The adopted
wind velocity distribution (β = 1) is a typical one for winds of massive stars, and more realistic thanin the simulation shown in Fig. 1. The ﬁgure shows the structure of colliding winds within r = a. The
upper and lower panels show the density and temperature distribution on the orbital plane, respectively.
From left to right, the orbital phases are 32 days prior to the periastron and 1 day and 30 days after
the periastron. In the ﬁgure, we see the eﬀect of radiative cooling (and that of winds signiﬁcantly
slower than the terminal velocities). Around the periastron, about half of the bow shocked region near
the primary cooled down to the ﬂoor temperature. At present, our code cannot reliably calculate theradiative cooling down to the temperature that allows dust formation. However, the result shown in
Fig. 2 may suggest the possibility that dust can form in the post shock region near periastron.
3.3.

Wolf-Rayet binary WR 140

WR 140 is an extensively studied, Wolf-Rayet binary consisting of a WC7 star and an O4-5V star. It
shows transient dust formation at the periastron. Both radial velocities and IR and X-ray light curves
vary with 7.9 yr periodicity. The X-ray light curve shows a deep minimum at the periastron, which
resemble the X-ray behavior of η Carinae.
The simulation parameters we adopted for this system are summarized in Table 1. We took thewind velocities and mass loss rates from model A of Zhekov & Skinner (2000), which makes the wind
momentum ﬂux ratio η  0.04. The value of η, however, seems to have a large ambiguity (0.04-0.2),
depending on the method used to derive it in previous studies. Below, as in the case of η Car, we ﬁrst
plot the large-scale density structure of the colliding winds of WR 140 and then show how (in)eﬃcient
the radiative cooling in this system is.
Figure 3 shows the density distribution of the colliding winds within r = 100a from WR 140. The
orbital phase of the plots is 200 days after the periastron passage. The format of the ﬁgure is the sameas that of Fig. 1(a). From Fig. 3, we note that in WR 140 it is the lower-density, primary (O4-5V)
wind that carves out a spiral cavity in the denser wind from the secondary (WC7). Because of thevery-high orbital eccentricity, interaction surface is strongly asymmetric, as of the winds of η Car: thecavities are very thin on the periastron side, whereas on the apastron side they occupy a large volume
separated by thin dense shells. A closer look, however, reveals diﬀerences caused by the diﬀerences in
the wind momentum ratio and the speed of the slower wind: the shock opening angle is wider and thespiral structure is more tightly wound in η Car than in WR 140. These diﬀerences are likely to aﬀect
the observational appearance of these systems.
Fig. 4 presents snapshots from a simulation where the radiative cooling is taken into account. It shows

147

Structure of Wind-Wind Collision Shocks in Massive Binaries

5

Fig. 4. Density (upper panels) and temperature (lower panels) distribution of wind-wind collision shocks in theorbital plane in the vicinity of the binary stars of WR 140 (r ≤ a). From left to right, the orbital phases are t = −30 d,
t = 0, and t = +30 d. Optically-thin radiative cooling is taken into account, with the floor temperature at 42,000 K.
Unperturbed wind has constant (β = 0) velocity distribution.

density distribution of colliding winds within r = a from WR 140. In this simulation, the β = 0 (i.e.,
constant) wind velocity distribution is adopted. The format of the ﬁgure is the same as that of Fig. 2.
From left to right, the orbital phases are 30 days prior to the periastron, the periastron, and 30 days
after it. In contrast to the case of η Car, we see little eﬀect of radiative cooling on the shocked regions.
The temperature of the apex shock is still as high as  107 K. This is because the primary’s wind speed
is much higher and the density is much lower than those in η Car. Thus, the origin of dust emission
from WR 140 is still an open question.
4.

Concluding Remarks

Constructing 3-D dynamical models is essential to understanding of the behavior of colliding-windbinaries. This is particularly the case for those with high orbital eccentricities. Thanks to the recent
development of 3-D dynamical modeling, we now have much better understanding of the hydrodynamic
interactions of wind-wind collision shocks than we did several years ago. Formation of dust in strong
shocks of colliding winds, however, is still a tough challenge to numerical simulations. Our results on
the supermassive binary η Carinae and the Wolf-Rayet binary WR 140 show that shocked regions have
too high temperature to allow the dust formation, although there is a hint that the temperature in thepost shock region can become low if the wind is very dense and slow as in η Car.

References
Bate M.R., Bonnell I.A., Price N.M., 1995, MNRAS, 285, 33
Canto, J., Raga, A. C., & Wilkin, F. P. 1996, ApJ, 469, 729
Corcoran, M. F., Hamaguchi, K., Pittard, J. M., Russell, C. M. P., Owocki, S. P., Parkin, E. R., Okazaki,
A. T., ApJ, 725, 1528
Damineli, A., 1996, ApJ, 460, L49
Ferland G.J., 1996, CLOUDY: 90.01
Gayley, K. G. 2009, ApJ, 703, 89
Lemaster, M. N., Stone, J. M., Gardiner, T. A., 2007, ApJ, 662, 582

148

6

A. T. Okazaki

Okazaki A.T., Bate M.R., Ogilvie G.I., Pringle J.E., 2002, MNRAS, 337, 967
Okazaki, A. T., Owocki, S. P., Russell, C. M. P., Corcoran, M. F., 2008, MNRAS: Letters, 388, L39
Parkin, E. R., Pittard, J. M., 2008, MNRAS, 388, 1047
Parkin, E. R., Pittard, J. M., Corcoran, M. F., Hamaguchi, K., 2011, ApJ, 726, 105
Pittard, J. M.; Stevens, I. R., 1999, In: IAUS 193, Wolf-Rayet Phenomena in Massive Stars and Starburst
Galaxies, ed. K. A. van der Hucht, G. Koenigsberger, P. R. J. Eenens (ASP, San Francisco), 386
Pittard, J. M., 2009, MNRAS, 396, 1743
Smith, N., 2010, MNRAS, 402, 145
Stevens, I. R., Blondin, J. M., & Pollock, A. 1992, ApJ, 386, 265
Tuthill, P. G., Monnier, J. D., Lawrance, N., Danchi, W. C., Owocki, S. P., Gayley, K. G., ApJ, 675,
698
Walder, R., Folini, D., 2003, In: IAUS 212, A Massive Star Odyssey: From Main Sequence to Supernova,
ed. K. van der Hucht, A. Herrero, E. César (ASP, San Francisco), 139
Zhekov S.A. & Skinner S.L., 2000, ApJ, 538, 808