We study GeV emission from gamma-ray binaries by assuming that the compact object is a young pulsar. We assume that the relativistic unshocked pulsar wind with Lorentz factor of 10^4~5 can produce the GeV emission by the inverse-Compton scattering process in the dense soft-photon field of the companion star. The travel distance of the unshocked pulsar wind that moves toward the observer depends on the orbital phase of the pulsar. We discuss that the orbital modulation of the GeV emission is a result of combination of the effects of the travel distance of the unshocked pulsar wind and of the anisotropic soft-photon field of the companion star. In this paper, we study how the effect of the travel distance of the unshocked pulsar wind affects to the orbital modulation of GeV emission. We apply our scenario to two gamma-ray binaries, LMC P3 and 4FGL J1405.5-6119. We find that with the suggested system parameters of LMC P3, the observed amplitude of the orbital modulation and the peak width are more consistent with the model light curve by taking into account the effect of the travel distance. For LMC P3, we analyze the GeV spectrum with 8-years Fermi-LAT data and discuss the broadband emission process in X-ray to TeV energy bands. We predict a possible system geometry for 4FGL J1405.5-6119 by fitting the GeV light curve. Key words: shock waves, gamma-rays: stars, stars: massive stars: mass-loss, stars: neutron

Substructures in PPDs, whose ubiquity was unveiled by recent ALMA observations, are widely discussed regarding their possible origins. We carry out global full magnetohydrodynamic (MHD) simulations in axisymmetry, coupled with self-consistent ray-tracing radiative transfer, thermochemistry, and non-ideal MHD diffusivities. The abundance profiles of grains are also calculated based on the global dust evolution calculation, including sintering effects. We found that dust size plays a crucial role in the ring formation around the snow lines of protoplanetary disks (PPDs) through the accretion process. Disk ionization structures and thus tensorial conductivities depend on the size of grains.When grains are significantly larger than PAHs, the non-ideal MHD conductivities change dramatically across each snow line of major volatiles, leading to a sudden change of the accretion process across the snow lines and the subsequent formation of gaseous rings/gaps there. On the other hand,the variations of conductivities are a lot less with only PAH sized grains in disks and then these disks retain smoother radial density profiles across snow lines.

Recent high resolution observations unveil ring structures in circumstellar disks. The origin of these rings has been widely investigated under various theoretical scenarios. In this work we perform global 3D non-ideal MHD simulations including effects from both Ohmic resistivity and ambipolar diffusion (AD) to model the HL Tau disk. The non-ideal MHD diffusion profiles are calculated based on the global dust evolution calculation including sintering effects. Disk ionization structure changes dramatically across the snow line due to the change of dust size distribution close to snow line of major volatiles. We find that accretion is mainly driven by disk wind. Gaps and rings can be quickly produced from different accretion rates across snow line. Furthermore, ambipolar diffusion (AD) leads to highly preferential accretion at midplane, followed by magnetic reconnection. This results a local zone of decretion that drains of mass in the field reconnection area, which leaves a gap and an adjacent ring just outside it. Overall, under the favorable condition, both snow lines and non-ideal MHD effects can lead to gaseous gaps and rings in protoplanetary disks.

We use the Fisher information matrix method to calculate the parameter estimation accuracy of inspiraling supermassive black holes binaries for TianQin, a space-borne laser interferometric detector aimed at detecting gravitational waves in the millihertz frequency band. The `restricted' post-Newtonian waveform in which third order post-Newtonian (3PN) phase including spin effects (spin-orbit $\beta$ and spin-spin $\sigma$) and first-order eccentricity contribution is employed. Monte Carlo simulations using $10^3$ binaries for mass pairs with component masses in the range of $({10^5},{10^7}){M_ \odot }$ and cosmological redshift $z=0.5$ show that the medians of the root-mean-square error distributions for the chirp mass $M_c$ and symmetric mass ratio $\eta$ are in the range of $\sim 0.02\% - 0.7\% $ and $\sim 4\% - 8\% $, respectively. The luminosity distance $D_L$ can be determined to be $\sim 1\% - 3\% $, and the angular resolution of source $\Delta \Omega $ is better than 12 deg$^2$. The corresponding results for $z=1.0$ and $2.0$, which are deteriorated with the decreasing of the signal-to-noise ratio, have also been given. We show that adding spin parameters degrades measurement accuracy of the mass parameters (${M_c}$, $\eta$), and the time and the orbital phase of coalescence ($t_c$, $\phi _c$); the inclusion of the first-order eccentricity correction to the phase worsens the estimation accuracy comparing with the circular cases. We also show the effects of post-Newtonian order on parameter estimation accuracy by comparing the results based on second order and third order post-Newtonian phases. Moreover, we calculate the horizon distance of supermassive black hole binaries for TianQin.

On 27 Apr 2015, when 67P/C-G was at 1.76 au from the Sun and moving towards perihelion, the OSIRIS and VIRTIS-M instruments on Rosetta observed the evolving dust and gas coma during a complete rotation of the comet. We aim to characterize the dust, H2O and CO2 gas spatial distribution in the inner coma. To do this we performed a quantitative analysis of the release of dust and gas and compared the observed H2O production rate with the one calculated using a thermo-physical model. For this study we selected OSIRIS WAC images at 612 nm (dust) and VIRTIS-M image cubes at 612 nm, 2700 nm (H2O) and 4200 nm (CO2). We measured the average signal in a circular annulus, to study spatial variation around the comet, and in a sector of the annulus, to study temporal variation in the sunward direction with comet rotation, both at a fixed distance of 3.1 km from the comet centre. The spatial correlation between dust and water, both coming from the sun-lit side of the comet, shows that water is the main driver of dust activity in this time period. The spatial distribution of CO2 is not correlated with water and dust. There is no strong temporal correlation between the dust brightness and water production rate as the comet rotates. The dust brightness shows a peak at 0deg sub-solar longitude, which is not pronounced in the water production. At the same epoch, there is also a maximum in CO2 production. An excess of measured water production, with respect to the value calculated using a simple thermo-physical model, is observed when the head lobe and regions of the Southern hemisphere with strong seasonal variations are illuminated. A drastic decrease in dust production, when the water production (both measured and from the model) displays a maximum, happens when typical Northern consolidated regions are illuminated and the Southern hemisphere regions with strong seasonal variations are instead in shadow.

TianQin is a space-based laser interferometric gravitational wave detector aimed at detecting gravitational waves at low frequencies (0.1 mHz -- 1 Hz). It is formed by three identical drag-free spacecrafts in an equilateral triangular constellation orbiting around the Earth. The distance between each pair of spacecrafts is approximately $1.7 \times 10^{5} ~\rm{km}$. The spacecrafts are interconnected by infrared laser beams forming up to three Michelson-type interferometers. The detailed mission design and the study of science objectives for the TianQin project depend crucially on the orbit and the response of the detector. In this paper, we provide the analytic expressions for the coordinates of the orbit for each spacecraft in the heliocentric-ecliptic coordinate system to the leading orders. This enables a sufficiently accurate study of science objectives and data analysis, and serves as a first step to further orbit design and optimization. We calculate the response of a single Michelson detector to plane gravitational waves in arbitrary waveform which is valid in the full range of the sensitive frequencies. It is then used to generate the more realistic sensitivity curve of TianQin. We apply this model on a reference white-dwarf binary as a proof of principle.

Systems with tightly-packed inner planets (STIPs) are very common. Chatterjee & Tan proposed Inside-Out Planet Formation (IOPF), an in situ formation theory, to explain these planets. IOPF involves sequential planet formation from pebble-rich rings that are fed from the outer disk and trapped at the pressure maximum associated with the dead zone inner boundary (DZIB). Planet masses are set by their ability to open a gap and cause the DZIB to retreat outwards. We present models for the disk density and temperature structures that are relevant to the conditions of IOPF. For a wide range of DZIB conditions, we evaluate the gap opening masses of planets in these disks that are expected to lead to truncation of pebble accretion onto the forming planet. We then consider the evolution of dust and pebbles in the disk, estimating that pebbles typically grow to sizes of a few cm during their radial drift from several tens of AU to the inner, $\lesssim1\:$AU-scale disk. A large fraction of the accretion flux of solids is expected to be in such pebbles. This allows us to estimate the timescales for individual planet formation and entire planetary system formation in the IOPF scenario. We find that to produce realistic STIPs within reasonable timescales similar to disk lifetimes requires disk accretion rates of $\sim10^{-9}\:M_\odot\:{\rm yr}^{-1}$ and relatively low viscosity conditions in the DZIB region, i.e., Shakura-Sunyaev parameter of $\alpha\sim10^{-4}$.

The processes that led to the formation of the planetary bodies in the Solar System are still not fully understood. Using the results obtained with the comprehensive suite of instruments on-board ESA's Rosetta mission, we present evidence that comet 67P/Churyumov-Gerasimenko likely formed through the gentle gravitational collapse of a bound clump of mm-sized dust aggregates ("pebbles"), intermixed with microscopic ice particles. This formation scenario leads to a cometary make-up that is simultaneously compatible with the global porosity, homogeneity, tensile strength, thermal inertia, vertical temperature profiles, sizes and porosities of emitted dust, and the steep increase in water-vapour production rate with decreasing heliocentric distance, measured by the instruments on-board the Rosetta spacecraft and the Philae lander. Our findings suggest that the pebbles observed to be abundant in protoplanetary discs around young stars provide the building material for comets and other minor bodies.

We present a statistical analysis of the distribution of large scale topographic features on comet 67P/Churyumov-Gerasimenko. We observe that the cumulative cliff height distribution across the surface follows a power law with a slope equal to -1.69 +- 0.02. When this distribution is studied independently for each region, we find a good correlation between the slope of the power law and the orbital erosion rate of the surface. For instance, the northern hemisphere topography is dominated by structures on the 100~m scale while the southern hemisphere topography, illuminated at perihelion, is dominated by 10~m scale terrain features. Our study suggest that the current size of a cliff is controlled not only by material cohesion but by the dominant erosional process in each region. This observation can be generalized to other comets, where we argue that primitive nuclei are characterized by the presence of large cliffs with a cumulative height power index equal to or above -1.5, while older, eroded cometary surfaces have a power index equal to or below -2.3. In effect, our model shows that a measure of the topography provides a quantitative assessment of a comet's erosional history, i.e. its evolutionary age.

Images of the nucleus and the coma (gas and dust) of comet 67P/Churyumov- Gerasimenko have been acquired by the OSIRIS (Optical, Spectroscopic, and Infrared Remote Imaging System) cameras since March 2014 using both the Wide Angle Camera (WAC) and the Narrow Angle Camera (NAC). We use images from the NAC camera to study a bright outburst observed in the southern hemisphere on July 29, 2015. The high spatial resolution of the NAC is needed to localize the source point of the outburst on the surface of the nucleus. The heliocentric distance is 1.25 au and the spacecraft-comet distance is 186 km. Aiming to better understand the physics that led to the outgassing, we used the Direct Simulation Monte Carlo (DSMC) method to study the gas flow close to the nucleus and the dust trajectories. The goal is to understand the mechanisms producing the outburst. We reproduce the opening angle of the outburst in the model and constrain the outgassing ratio between the outburst source and the local region. The outburst is in fact a combination of both gas and dust, in which the active surface is approximately 10 times more active than the average rate found in the surrounding areas. We need a number of dust particles 7.83 $\times$ 10$^{11}$ - 6.90 $\times$ 10$^{15}$ (radius 1.97 - 185 {\mu}m), which corresponds to a mass of dust 220 - 21 $\times$ 10$^{3}$kg.

Beginning in March 2014, the OSIRIS (Optical, Spectroscopic, and Infrared Remote Imaging System) cameras began capturing images of the nucleus and coma (gas and dust) of comet 67P/Churyumov-Gerasimenko using both the wide angle camera (WAC) and the narrow angle camera (NAC). The many observations taken since July of 2014 have been used to study the morphology, location, and temporal variation of the comet's dust jets. We analyzed the dust monitoring observations shortly after the southern vernal equinox on May 30 and 31, 2015 with the WAC at the heliocentric distance Rh = 1.53 AU, where it is possible to observe that the jet rotates with the nucleus. We found that the decline of brightness as a function of the distance of the jet is much steeper than the background coma, which is a first indication of sublimation. We adapted a model of sublimation of icy aggregates and studied the effect as a function of the physical properties of the aggregates (composition and size). The major finding of this article was that through the sublimation of the aggregates of dirty grains (radius a between 5 microm and 50 microm) we were able to completely reproduce the radial brightness profile of a jet beyond 4 km from the nucleus. To reproduce the data we needed to inject a number of aggregates between 8.5 x $10^{13}$ and 8.5 x $10^{10}$ for a = 5 microm and 50 microm respectively, or an initial mass of $H_2O$ ice around 22kg.

A new air-shower core-detector array (YAC : Yangbajing Air-shower Core-detector array) has been developed to measure the primary cosmic-ray composition at the "knee" energies in Tibet, China, focusing mainly on the light components. The prototype experiment (YAC-I) consisting of 16 detectors has been constructed and operated at Yangbajing (4300 m a.s.l.) in Tibet since May 2009. YAC-I is installed in the Tibet-III AS array and operates together. In this paper, we performed a Monte Carlo simulation to check the sensitivity of YAC-I+Tibet-III array to the cosmic-ray light component of cosmic rays around the knee energies, taking account of the observation conditions of actual YAC-I+Tibet-III array. The selection of light component from others was made by use of an artificial neural network (ANN). The simulation shows that the light-component spectrum estimated by our methods can well reproduce the input ones within 10\% error, and there will be about 30\% systematic errors mostly induced by the primary and interaction models used. It is found that the full-scale YAC and the Tibet-III array is powerful to study the cosmic-ray composition, in particular, to obtain the energy spectra of protons and helium nuclei around the knee energies.

Dust jets, i.e. fuzzy collimated streams of cometary material arising from the nucleus, have been observed in-situ on all comets since the Giotto mission flew by comet 1P/Halley in 1986. Yet their formation mechanism remains unknown. Several solutions have been proposed, from localized physical mechanisms on the surface/sub-surface (see review in Belton (2010)) to purely dynamical processes involving the focusing of gas flows by the local topography (Crifo et al. 2002). While the latter seems to be responsible for the larger features, high resolution imagery has shown that broad streams are composed of many smaller features (a few meters wide) that connect directly to the nucleus surface. We monitored these jets at high resolution and over several months to understand what are the physical processes driving their formation, and how this affects the surface. Using many images of the same areas with different viewing angles, we performed a 3-dimensional reconstruction of collimated jets, and linked them precisely to their sources on the nucleus. Results.We show here observational evidence that the Northern hemisphere jets of comet 67P arise from areas with sharp topographic changes and describe the physical processes involved. We propose a model in which active cliffs are the main source of jet-like features, and therefore the regions eroding the fastest on comets. We suggest that this is a common mechanism taking place on all comets.

The Kepler mission has discovered more than 4000 exoplanet candidates. Many are in systems with tightly packed inner planets. Inside-Out Planet Formation (IOPF) has been proposed to explain these systems. It involves sequential in situ planet formation at the local pressure maximum of a retreating dead zone inner boundary (DZIB). Pebbles accumulate at this pressure trap, which builds up a ring, and then a planet. The planet is expected to grow until it opens a gap, which helps to both truncate pebble accretion and induce DZIB retreat that sets the location of formation of the next planet. This simple scenario may be modified if the planet migrates significantly from its formation location. Thus planet-disk interactions play a crucial role in the IOPF scenario. We present numerical simulations that first assess migration of planets of various masses that are forming at the DZIB of an active accretion disk, where the effective viscosity rapidly increases in the radially inward direction. We find that the disk's torques on the planet tend to trap the planet at a location very close to the initial pressure maximum where it formed. We then study gap opening by these planets to assess at what mass a significant gap is created. Finally we present a simple model for DZIB retreat due to penetration of X-rays from the star to the disk midplane. Overall, these simulations help to quantify both the mass scale of first,"Vulcan," planet formation and the orbital separation to the location of second planet formation.

The Kepler-discovered Systems with Tightly-packed Inner Planets (STIPs), typically with several planets of Earth to super-Earth masses on well-aligned, sub-AU orbits may host the most common type of planets, including habitable planets, in the Galaxy. They pose a great challenge for planet formation theories, which fall into two broad classes: (1) formation further out followed by inward migration; (2) formation in situ, in the very inner regions of the protoplanetary disk. We review the pros and cons of these classes, before focusing on a new theory of sequential in situ formation from the inside-out via creation of successive gravitationally unstable rings fed from a continuous stream of small (~cm-m size) "pebbles," drifting inward via gas drag. Pebbles first collect at the pressure trap associated with the transition from a magnetorotational instability (MRI)-inactive ("dead zone") region to an inner, MRI-active zone. A pebble ring builds up that begins to dominate the local mass surface density of the disk and spawns a planet. The planet continues to grow, most likely by pebble accretion, until it becomes massive enough to isolate itself from the accretion flow via gap opening. This reduces the local gas density near the planet, leading to enhanced ionization and a retreat of the dead zone inner boundary. The process repeats with a new pebble ring gathering at the new pressure maximum associated with this boundary. We discuss the theory's predictions for planetary masses, relative mass scalings with orbital radius, and minimum orbital separations, and their comparison with observed systems. Finally, we discuss open questions, including potential causes of diversity of planetary system architectures, i.e., STIPs versus Solar System analogs.

Aiming at the observation of cosmic-ray chemical composition at the "knee" energy region, we have been developinga new type air-shower core detector (YAC, Yangbajing Air shower Core detector array) to be set up at Yangbajing (90.522$^\circ$ E, 30.102$^\circ$ N, 4300 m above sea level, atmospheric depth: 606 g/m$^2$) in Tibet, China. YAC works together with the Tibet air-shower array (Tibet-III) and an underground water cherenkov muon detector array (MD) as a hybrid experiment. Each YAC detector unit consists of lead plates of 3.5 cm thick and a scintillation counter which detects the burst size induced by high energy particles in the air-shower cores. The burst size can be measured from 1 MIP (Minimum Ionization Particle) to $10^{6}$ MIPs. The first phase of this experiment, named "YAC-I", consists of 16 YAC detectors each having the size 40 cm $\times$ 50 cm and distributing in a grid with an effective area of 10 m$^{2}$. YAC-I is used to check hadronic interaction models. The second phase of the experiment, called "YAC-II", consists of 124 YAC detectors with coverage about 500 m$^2$. The inner 100 detectors of 80 cm $\times $ 50 cm each are deployed in a 10 $\times$ 10 matrix from with a 1.9 m separation and the outer 24 detectors of 100 cm $\times$ 50 cm each are distributed around them to reject non-core events whose shower cores are far from the YAC-II array. YAC-II is used to study the primary cosmic-ray composition, in particular, to obtain the energy spectra of proton, helium and iron nuclei between 5$\times$$10^{13}$ eV and $10^{16}$ eV covering the "knee" and also being connected with direct observations at energies around 100 TeV. We present the design and performance of YAC-II in this paper.

Static and dynamic properties of low density outer envelopes of neutron stars are calculated within the nonlinear magnetic Thomas-Fermi model, assuming degenerate electrons. A novel domain decomposition enables proper description of lattice symmetry and may be seen as a prototype for the general class of problems involving nonlinear charge screening of periodic, quasi-low-dimensionality structures, e.g. liquid crystals. We describe a scalable implementation of the method using Hypre. Phase velocity of long wavelength transverse phonons is found to be a factor of 5-7 larger than in the corresponding Coulomb crystal model, which could have implications for low temperature phonon-mediated thermal conductivity. Other findings include $c'<0$ elastic instabilities for both bcc and fcc lattices, reminiscent of the situation in some light actinides, and suggestive of a symmetry-lowering transition to a tetragonal or orthorhombic lattice.

Inside-Out Planet Formation (IOPF; Chatterjee & Tan 2014, hereafter CT14) is a scenario for sequential in situ planet formation at the pressure traps of retreating dead zone inner boundaries (DZIBs) motivated to explain the many systems with tightly packed inner planets (STIPs) discovered by Kepler. The scenario involves build-up of a pebble-dominated protoplanetary ring, supplied by radial drift of pebbles from the outer disk. It may also involve further build-up of planetary masses to gap-opening scales via continued pebble accretion. Here we study radial drift & growth of pebbles delivered to the DZIB in fiducial IOPF disk models.

Rate constants for the N + C2 reaction have been measured in a continuous supersonic flow reactor over the range 57 K to 296 K by the relative rate technique employing the N + OH - H + NO reaction as a reference. Excess concentrations of atomic nitrogen were produced by the microwave discharge method and C2 and OH radicals were created by the in-situ pulsed laser photolysis of precursor molecules C2Br4 and H2O2 respectively. In parallel, quantum dynamics calculations were performed based on an accurate global potential energy surfaces for the three lowest lying quartet states of the C2N molecule. The 14A" potential energy surface is barrierless, having two deep potential wells corresponding to the NCC and CNC intermediates. Both the experimental and theoretical work show that the rate constant decreases to low temperature, although the experimentally measured values fall more rapidly than the theoretical ones except at the lowest temperatures. Astrochemical simulations indicate that this reaction could be the dominant source of CN in dense interstellar clouds.

The China Dark matter Experiment collaboration reports the first experimental limit on WIMP dark matter from 14.6 kg-day of data taken with a 994 g p-type point-contact germanium detector at the China Jinping underground Laboratory where the rock overburden is more than 2400 m. The energy threshold achieved was 400 eVee. According to the 14.6 kg-day live data, we placed the limit of N= 1.75 * 10^{-40} cm^{2} at 90% confidence level on the spin-independent cross-section at WIMP mass of 7 GeV before differentiating bulk signals from the surface backgrounds.

Rate constants for the potentially important interstellar N(4S) + CH(X2{\Pi}r) reaction have been measured in a continuous supersonic flow reactor over the range 56 K < T < 296 K using the relative rate technique employing both the N(4S) + OH(X2{\Pi}i) and N(4S) + CN(X2{\Sigma}+) reactions as references. Excess concentrations of atomic nitrogen were produced by the microwave discharge method upstream of the Laval nozzle and CH and OH radicals were created by the in-situ pulsed laser photolysis of suitable precursor molecules. In parallel, quantum dynamics calculations of the title reaction have been performed based on accurate global potential energy surfaces for the 13A' and 13A" states of HCN and HNC, brought about through a hierarchical construction scheme. Both adiabatic potential energy surfaces are barrierless, each one having two deep potential wells suggesting that this reaction is dominated by a complex-forming mechanism. The experimental and theoretical work are inexcellent agreement, predicting a positive temperature dependence of the rate constant, in contrast to earlier experimental work at low temperature. The effects of the new low temperature rate constants on interstellar N2 formation are tested using a dense cloud model, yielding N2 abundances 10-20 % lower than previously predicted.

The CDEX Collaboration has been established for direct detection of light dark matter particles, using ultra-low energy threshold p-type point-contact germanium detectors, in China JinPing underground Laboratory (CJPL). The first 1 kg point-contact germanium detector with a sub-keV energy threshold has been tested in a passive shielding system located in CJPL. The outputs from both the point-contact p+ electrode and the outside n+ electrode make it possible to scan the lower energy range of less than 1 keV and at the same time to detect the higher energy range up to 3 MeV. The outputs from both p+ and n+ electrode may also provide a more powerful method for signal discrimination for dark matter experiment. Some key parameters, including energy resolution, dead time, decay times of internal X-rays, and system stability, have been tested and measured. The results show that the 1 kg point-contact germanium detector, together with its shielding system and electronics, can run smoothly with good performances. This detector system will be deployed for dark matter search experiments.

Weakly Interacting Massive Particles (WIMPs) are the candidates of dark matter in our universe. Up to now any direct interaction of WIMP with nuclei has not been observed yet. The exclusion limits of the spin-independent cross section of WIMP-nucleon which have been experimentally obtained is about 10^{-7}pb at high mass region and only 10^{-5}pb} at low mass region. China Jin-Ping underground laboratory CJPL is the deepest underground lab in the world and provides a very promising environment for direct observation of dark matter. The China Dark Matter Experiment (CDEX) experiment is going to directly detect the WIMP flux with high sensitivity in the low mass region. Both CJPL and CDEX have achieved a remarkable progress in recent two years. The CDEX employs a point-contact germanium semi-conductor detector PCGe whose detection threshold is less than 300 eV. We report the measurement results of Muon flux, monitoring of radioactivity and Radon concentration carried out in CJPL, as well describe the structure and performance of the 1 kg PCGe detector CDEX-1 and 10kg detector array CDEX-10 including the detectors, electronics, shielding and cooling systems. Finally we discuss the physics goals of the CDEX-1, CDEX-10 and the future CDEX-1T detectors.

Unparticle $\U$ with scaling dimension $d_\U$ has peculiar thermal properties due to its unique phase space structure. We find that the equation of state parameter $\omega_\U$, the ratio of pressure to energy density, is given by $1/(2d_\U +1)$ providing a new form of energy in our universe. In an expanding universe, the unparticle energy density $\rho_\U(T)$ evolves dramatically differently from that for photons. For $d_\U >1$, even if $\rho_\U(T_D)$ at a high decoupling temperature $T_D$ is very small, it is possible to have a large relic density $\rho_\U(T^0_\gamma)$ at present photon temperature $T^0_\gamma$, large enough to play the role of dark matter. We calculate $T_D$ and $\rho_\U(T^0_\gamma)$ using photon-unparticle interactions for illustration.