\section{Introduction}

Large interferometer telescopes such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Very Large Telescope Interferometer (VLTI) have recently facilitated an increasingly detailed examination of circumstellar disks around young stars. They have the ability to detect molecular emission spectra, as well as dust continuum emission, at very high angular resolutions even at the millimeter and sub-millimeter wavelengths. The much higher angular resolution, when compared to a single telescope at similar wavelengths, allows us to probe closer to the central star, and study the spatial relationships between gas and dust on a finer scale than ever before. 

PDS 70 is a pre-main-sequence star with spectral type K5 and assumed to be part of the Centaurus association, estimated to be at a distance of 140 pc \citep{riaud06}. It has previously been imaged in both $H$-band polarized light using Subaru \citep{hashimoto12}, as well as more recently at 1.3 mm with the Sub Millimeter Array (SMA) \citep{hashimoto15}. In this paper we present a detailed examination of the circumstellar disk of PDS 70 in the 0.87 mm continuum, CO J=3--2 emission, and HCO\textsuperscript{+} J=4--3 emission based on ALMA observations. These data are at a much higher resolution than the SMA data, $0\farcs19 \times 0\farcs15$ compared to $0\farcs88 \times 0\farcs43$, and provide us with a good probe of the larger grains within the disk. They also allow us to examine the gas distribution and detect fine structures within disk gas. This is of particular interest to our study due to the existence of features which appear in the gas disk but are not present in the dust disk. PDS 70 has proven to be a system with a different structure in each of our continuum, CO J=3--2,  HCO\textsuperscript{+} J=4--3 measurements presented below. In addition to our observational data we present a Monte Carlo Radiative Transfer (MCRT) model of the disk using HOCHUNK3D \citep{whitney13} to simultaneously fit existing imagery as well as the spectral energy distribution (SED) of PDS 70. Following the discussion of the data and model, we also present additional SpeX and REM observations of PDS 70 as well as analysis of these data sets in the Appendix including gas accretion, mass, and age estimates.

\section{Observations and Data Reduction}

Observations were conducted on August 14 and 18, 2016 in the ALMA cycle 3 program (No. 2015.1.000888.S). We used 38 to 39 of the 12 m diameter antennas in the observation. The array configuration provides a maximum and minimum baseline length of 1.5 km and 15.1 m, respectively. Sky conditions were relatively stable with the precipitable water vapor between 0.45 mm and 0.82 mm. The system temperature was between 90 K and 280 K during observations. The correlator was tuned at the center frequency of 345.761 GHz for CO J=3–-2, 356.698 GHz for HCO\textsuperscript{+} J=4–-3, and 343.955 GHz and 355.497 GHz for continuum observations in dual polarization mode. The effective bandwidths for both emission lines are 468.75 MHz and their frequency resolution is 244.141 kHz, corresponding to approximately 0.21 km s\textsuperscript{-1} in the velocity resolution.

Data reduction was performed using version 4.7 of the Common Astronomy Software Applications package (CASA). A quasar J1427-4206 was used for bandpass, flux, and phase calibrations, and another quasar J1407-4302 was selected as a second calibrator. We used the task CLEAN in CASA and a continuum image was generated by combining the line-free channels in all of the spectral windows. The self-calibration was performed by applying the obtained continuum image as a model, and gain calibration was repeated until the rms reached minimum. The gain table obtained after the self-calibration of the continuum data was applied to both the CO J=3--2 and HCO\textsuperscript{+} J=4--3 to generate self-calibrated visibilities. Because the line and continuum observations were simultaneously conducted, and the atmospheric variation equally affects the phase variation of line and continuum data, the calibration table obtained from the continuum data can be applied to line emission calibration. Natural weighting was applied in both continuum and line imaging for better sensitivity. The achieved rms was 0.0551 mJy beam$^{-1}$ for the continuum, 1.204 mJy beam$^{-1}$ for CO J=3--2, and 1.118 mJy beam$^{-1}$ for HCO\textsuperscript{+} J=4--3. The field of view was approximately 18$\arcsec$ and the final synthesized beam of approximately $0\farcs19\times0\farcs15$ for the continuum image and CO J=3--2 map, and $0\farcs18\times0\farcs15$ for the HCO\textsuperscript{+} J=4--3 map after calibration and flagging as denoted in Figures~\ref{fig:continuum},~\ref{fig:gas_data0}. The total integration time on source after flagging is about 56 minutes. 


\section{Results}

\subsection{0.87 mm Continuum Measurements}

The 0.87 mm continuum image gives us a good measure of the radial extent of the dust disk, the flux density, and reveals some other interesting features. The major axis of the disk lies at $\sim$159$\degr$ East of North, measured via ellipse fitting of the brightest emission along the ring of the disk, (Figure~\ref{fig:continuum}) and measures a radial distance of $0\farcs42 \pm 0\farcs05$ ($\sim$60 au) from the center at the edge and $0\farcs79 \pm 0\farcs05$ ($\sim$110 au) at the outer edge. Both edges were found by measuring a distance equal to half the width of the beam size from the inner and outer edges of the disk that are detectable above noise. The measured flux density is 233 $\pm$ 3 mJy integrated at a radius of $1\farcs00$, past which there is no detectable emission above $\sim 2 \sigma$. In Figure~\ref{fig:continuum} the disk brightness is clearly not symmetrical about the minor axis, similar to what is seen in \citet{kraus17}, with the NW side $\sim$15$\%$ brighter than the SE side.

There also appears to be material visible near the center of the disk in the continuum, seen more clearly in the log scaled image (Figure~\ref{fig:continuum},  left), which is separated from the outer disk by a radial gap. It is slightly larger than the resolution of the measurement, extending out approximately 0.2$\arcsec$ where the resolution of the instrument is $0\farcs19 \times 0\farcs15$ and is likely an inner disk component.

\subsubsection{HCO\textsuperscript{+} J=4--3 Measurements}

The HCO\textsuperscript{+} J=4--3 moment 0 emission (Figure~\ref{fig:gas_data0}) reveals a similar gap to what is seen in the dust continuum (Figure~\ref{fig:continuum}), though it extends to only $\sim33$ au while the outer edge of the HCO\textsuperscript{+} J=4--3 moment 0 emission extends to $\sim160$ au. The gas disk also extends out farther than the dust though not as far as in many very young pre-main-sequence (PMS) stars such as DM Tau \citep{oberg10}. It contains a bright blob ($\sim$ 10 $\sigma$) to the N, as well as a bridge structure connecting the inner region and the outer disk in the SW, and a gap due W which lies directly adjacent to the bridge. There is also evidence of an azimuthal gap or depleted region to the NE.

The bridge structure is likely to be real, with a signal-to-noise ratio of $\sim$ 5. In order to help verify this we generated a rough model of the bridge and azimuthal gap (Figure~\ref{fig:Simulation}). This was produced using the CASA 5.2 task simobserve with the thermal noise set to default and the antennas set to the ALMA cycle 4 configuration, alma.cycle4.3.cfg. The seed image for the model is a simple, two-component system generated using HOCHUNK3D with the bridge and depleted region added to the seed using image editing software. We found that in order to produce any asymmetric features to the degree of those seen in Figure~\ref{fig:gas_data0}, the inclusion of dimmer region to reproduce the gap and a bridge connecting the inner and outer disks was needed.

The bridge feature itself appears to be mass transfer between the inner and outer disks similar to what is suggested for HD 142527 \citep{casassus13}. If this were the case we would expect to see a radial velocity component roughly perpendicular to the Keplerian motion of the disk which is not seen here (Figure~\ref{fig:gas_data12}). It has also been suggested that bridges such as this can be caused by a polar outflow of gas \citep{alves17}. Again however, this is likely not the case however as there appears to be no significant change in velocity in the region of the bridge, and it truncates at the outer disk instead of continuing to expel matter outside of the disk. Its close proximity to the western gap suggests that the two are related and could be due to large body formation which has been shown to cause both depleted regions as well as spiral arms in other disks \citep{edgar08, dong15, matter16}. Interestingly enough however, there is no evidence of the bridge feature or an azimuthal gap in the 0.87 mm continuum ring or in $H$-band polarized imagery \citep{dong12}. 

As can be seen in Figure~\ref{fig:gas_data0}, the inner HCO\textsuperscript{+} J=4--3 disk appears slightly offset from the inner dust disk. To verify this we created radial traces of the continuum and moment 0 gas data, shown in Figure~\ref{fig:sb_slices}. There also may be some asymmetry in the outer gas disk as the HCO\textsuperscript{+} J=4--3 moment 0 peak intensities lie at different distances from the center of the disk. In the HCO\textsuperscript{+} J=4--3 moment 1 and 2 maps, the central region was clipped out due to excess noise in the region (Figure~\ref{fig:clipping}). 

\subsubsection{CO J=3--2 Measurements}

The CO J=3--2 moment 0 emission lies interior to the continuum and HCO\textsuperscript{+} J=4--3 emission, exhibiting little to no gap between the inner and outer disk. Similar to the HCO\textsuperscript{+} J=4--3 moment 0 data this disk extends further out than the dust disk, to $\sim210$ au, but again not as much as much younger PMS stars. This is promising for planet formation similar to what is seen in our solar system as there is still mass within 40 au of the star. The emission contains a blob in the S similar to what is seen in HCO\textsuperscript{+} J=4--3, but almost directly opposite its location in the disk. 

Again, similar to what is seen in HCO\textsuperscript{+} J=4--3, there is a bar feature towards the center of the image. This, however, appears to extend to either side of the disk and, as CO J=3--2 is generally considered optically thick in protoplanetary disks, could be causing the dimmed regions seen at either end of the bridge. It appears the bridge feature lies within the plane of the disk as the moment 1 data shows no deviation from Keplerian rotation in the region of the bridge (Figure~\ref{fig:gas_data12}).

The CO moment 2 map contains a butterfly pattern surrounding the central region. This pattern exists in the same locations as the four bright regions in the moment 0 data, but are due to crowded isovelocity curves in the moment 1 data and are likely not related.


\subsection{Model}

\subsubsection{Initial Model}
For our modeling we used the MCRT code HOCHUNK3D \citep{whitney13}. This code allows for two distinct grain populations as well as an envelope which can be manipulated into a third distinct grain population if necessary. The locations of these populations can also be changed to produce radial gaps in the disk.

We used a model similar to those of \citet{hashimoto15} and \citet{dong12} as our starting point. For the star we used a pre-main-sequence star with spectral type K5, a radius of 1.39 R$\astrosun$, and a temperature of 4400 K as suggested by \citet{gregorio02}. We also use a mass of 0.82 M$\astrosun$, estimated by \citet{riaud06} using spectral type, and further supported by dynamical fitting \citep{hashimoto15}. We then constructed a disk consisting of two rings with a large gap in the middle, a mass of 4.5 M\textsubscript{J}, an inclination of 50$\degr$ \citep{hashimoto12}, and a large grain to total grain mass ratio of 0.9667 in order to decrease the optical thickness of the small grains and increase the 10 and 20 $\mu$m silicate peaks. We also used the same grain types as the \citet{hashimoto15} model for the large and small grains which make up the disk.

Our initial model consisted of one that included both grain types in the inner disk which extended from the sublimation radius $\sim$0.05 au to 1 au which was used to reproduce the near IR emission seen in the SED. This was followed by a gap which extended to 65 au for the smaller grains and 80 au for the larger grains based on the $H$-band data \citep{dong12}, the 1.3 mm data \citep{hashimoto15}, and a distance of 140 pc \citep{riaud06}. This reproduced a fit to both the SED and 1.3 mm image seen in \citet{hashimoto15} as well as \citet{dong12}. 

\subsubsection{Model Fitting} 

Since the population of small grains is constrained by the character of the silicate emission bands near the 10 and 20 $\mu m$ wavelengths, we also include data obtained with the $Spitzer$ InfraRed Spectrograph (IRS) in our analysis. Both of these features are the result of warm silicate emission close to the star and required modifying the initial model to provide a reasonable fit. 

It is important to note also that the AKARI 18 $\micron$ point disagrees with the $Spitzer$ spectrum despite the $13.9-23.6$ $\micron$ band pass of the filter, and could indicate that either there is some variability of the disk or one of the data sets is incorrect. The grains which produce this peak lie approximately 0.5-3 au from the star which, assuming that the star is 0.82 M$\astrosun$, suggests orbital periods $\sim$0.39$-$5.7 years. It is possible therefore, that the disk could have changed between the two data sets. Additionally, it has been suggested that thermal instabilities due to active accretion could result in variability \citep{sitko12}. However PDS 70 seems to exhibit no accretion (Appendix A).

A second difference is the inclusion of the newly-acquired of ALMA 0.87 mm continuum data described in Section 2. The ALMA data have a much higher resolution than the SMA data used by \citet{hashimoto15}, which allows for much better spatial estimates of the disk. These data exhibit a different, $\sim$60 au, gap radius and lack the highly asymmetric spatial distributions in the SMA dust emission. This is in agreement with the gap radius found in scattered light \citep{dong12} which suggests that the apparent difference in gap radius for the two grain populations is not real

To fit the model with these new constraints we included a less dense, though still optically thick, inner disk disk comprised of only small grains which extends from the sublimation radius to 15 au. We also raised the disk mass from 4.5 M\textsubscript{J} to $\sim7.3$ M\textsubscript{J} to fit the ALMA 0.87 mm point and raised the large grains to total grains ratio to 0.995.  This produced very strong silicate peaks, but does not reach the depth of the dips before the 10 and 20 $\micron$ peaks. The inner disk is followed by a gap extending to 60 au and an outer disk extending from 60 au to 110 au, based on the 0.87 mm continuum image. The higher spatial resolution of the ALMA data also allows for better comparison of the disk inclination to the model. We found the apparent inclination of the disk is best reproduced with a model inclination of 45$\degr$ from face on.

In order to probe the optical depth of the outer disk we examined the thickness of the disk, defined by $z = C_1 r^b$ where z is the scale height (thickness) of the disk, $C_1$ is a constant, r is the radial distance from the star, and b is the flaring exponent \citep{whitney03a}. Similarly the radial density is defined by $\rho = C_2 r^{-a}$ where $\rho$ is the density of the disk, $C_2$ is a constant, r is the radial distance from the star, and a is the density exponent \citep{whitney03a}. To match the SED we find b = 1.15 $\pm$ 0.02 and a = 2.25 $\pm$ 0.02 for the large grains in our disk and b = 1.15 $\pm$ 0.02 and a = 1.5 $\pm$ 0.02 for the small grains. Based on these parameters we find that the optical depth of the outer disk is $\tau \sim100$ au which is larger than the thickness of the outer disk suggesting that the disk is mostly optically thin at 0.87 mm. These two exponents can also be used to calculate the surface density of the disk via the method described in \citet{whitney13} eq. 15 ($\Sigma = \Sigma_0 \varpi^{-p}$) where $\Sigma_0$ is a constant set by the mass of the disk, $\varpi$ is the cylindrical radius, and $p=\alpha-\beta$. The surface density profile for the disk can be seen in Figure~\ref{fig:sigma}.

Because of the higher resolution of our data and the differences seen between this and the SMA data \citep{hashimoto15}, such as the outer disk extent, we used the CO J=3--2 moment 1 data to find the dynamical mass of PDS 70. Here we assumed a disk inclination of 45$\degr$, based on our modeling, and similar to what is found in \citet{hashimoto15} we see that a mass of 0.6 - 0.8 M$\astrosun$ (Figure~\ref{fig:dynamical_mass}) produces a good fit to the data and is in agreement with the \citet{riaud06} estimation. This falls within the range of 0.5 - 1.1 M$\astrosun$ found from comparing the location of PDS 70 to isochrones (Appendix) derived in \citet{tognelli11} and seen in Figure~\ref{fig:TTS_mass}.

\subsubsection{Final Model}

To test the accuracy of our model, we compared model images to both the 0.87 mm continuum data (Figure~\ref{fig:model_87}) and the $H$-band data \citep{hashimoto12} (Figure~\ref{fig:model_hband}) after convolving with beam shape and point-spread function of these data respectively. In the case of the 0.87 mm model image we see a ratio of emission between the inner and outer disks which is in agreement with the ALMA image, \(\frac{Finner}{Fouter}\) $\sim0.004$ in both cases, which suggests the inner disk is essentially devoid of large grains. The $H$-band data reproduces the apparent dip in brightness seen in the NE side of the disk (Figure~\ref{fig:model_hband}). This dip is the result of co-planar shadowing from the the inner disk. The final model parameters can bee seen in Table~\ref{tbl:parameters}.

\section{Discussion}

\subsection{Asymmetry in 0.87 mm Continuum Image}
If we assume the 0.87 mm continuum is optically thin, section 3.2.2, the difference in brightness should correspond to a difference in mass density between the two locations, suggesting that there is more mass density in the NW side of the outer disk. There exists the possibility however, that the difference in brightness could be a result of apocenter glow \citep{pan16}, though we do not detect a pericenter offset in this system. Given the resolution of the image however, there could be an offset up to $\sim5$ au which would result in an apocenter glow of $\sim6.7\%$  according to the calculation proposed in \citet{pan16}.

\subsection{Differences in Gas and Dust Imagery}

Most of the gas features seen in the HCO\textsuperscript{+} J=4--3 and CO J=3--2 maps do not have an analogue in the 0.87 mm continuum or $H$-band imagery. The only notable exceptions are the excess flux in the NW of the continuum image, an inner gap, and an inner disk which is present in all available data.

\subsubsection{Azimuthal Gap in HCO\textsuperscript{+} J=4--3}

While some of these features are easily explained by known properties, such as the diminished, central region in the CO J=3--2 emission and the larger radial extent of both HCO\textsuperscript{+} J=4--3 and CO J=3--2 emission, of the gas the lack of an azimuthal gap in the continuum or $H$-band data is not. It is likely not a gravitational depletion of the disk by a compact body as there should be no difference in how the dust and gas are affected gravitationally. This suggests a more exotic cause for the apparent difference.

Options we have considered are a large stellar flare from a coronal mass ejection, an embedded planet, or shadowing from an inner optically thick blob of material. A large stellar flare could disturb the gas in a localized region of the disk and possibly move the dust in that region as well. This could neatly link the azimuthal gap seen in the HCO\textsuperscript{+} J=4--3 data as well as the large blob of mass seen in the 0.87 mm continuum.  A similar effect was seen in AU Mic where it is suggested that a stellar flare caused fast moving features in the dust of the disk \citep{boccaletti15}. However, though this could cause feature propagation in the dust disk, we would expect the density wave to travel around the disk or perhaps cause dust blowout from the impact. The difference in density from either scenario is not seen in the 0.87 mm continuum and there seems to be no difference in the location of the bright blob in the 1.3 mm continuum \citep{hashimoto15} and the 0.87 mm continuum.

An embedded planet could cause the local region around it to undergo various chemical changes due to heating from the planet \citep{cleeves15}. Though the heating would initially cause any frozen out gasses to sublimate, it may also catalyze complicated chemical interactions which could locally deplete the HCO resulting in an apparent gap in the gas disk. It has been shown that a Jupiter-sized planet can cause an effect in a radius of several au surrounding the planet \citep{cleeves15} which grows as the planet moves further from the star. Unfortunately a planet of this size would not be visible in the currently available data.

Shadowing from an inner disk would result in an apparent gap in the disk dust \citep{long17,stolker16} and should similarly effect the gas if the inner disk is optically thick enough. In this instance, however, the shadowing body would necessarily have an asymmetric distribution as there is no equivalent azimuthal gap on the western side of the HCO\textsuperscript{+} J=4--3 disk, only minor dimming. This could explain the apparent inner disk offset seen in the HCO\textsuperscript{+} J=4--3 moment zero data as it lies on the NE side of the star closely aligned with the apparent gap. If this is the source of the apparent gap, then we would expect the shadow to move as the shadowing body moves around the star. This should happen within several months given that the shadowing body would appear to be at a radius of $\sim5$ au. 

\subsubsection{Gas-poor inner disk region}

As previously mentioned the SED suggests the existence of dust close to sublimation radius of the star which is further supported by our modeling. In contrast, we see in our SpeX data (discussed more thoroughly in the Appendix) that the innermost regions of the disk, less than 1 au, seem to be devoid of gas as there is no net emission for hydrogen or CO gas (Figure~\ref{fig:LineExtractions}) above photospheric. There is however some Ca II emission which suggests that there is chromoshperic activity similar to what is seen in HD 163296 \citep{sitko08}. The lack of gas emission close to the star suggests that there is also no active accretion of gas onto the star.



\subsubsection{Bar across CO}

The CO J=3--2 moment 1 data reveals what appears to be a bar along the minor axis of the disk.  At opposite sides of the disk, directly adjacent to the bridge feature there appear to be dips in brightness in the CO moment 1 map. If the bridge feature is optically thick, as expected from CO gas, it could in turn shadow the gas in the outer disk. This however raises the question as to how the bridge structure could have formed. It could be a mass transfer between the inner and outer disk, though over time we would expect that the bridge would not remain solely along the minor axis, but instead would smear out due to Keplerian motion. We would also expect that some dust would be transferred as well which is not seen in the continuum data. This suggests it could lie in a plane nearly perpendicular the plane of the outer disk. However, as mentioned in Section 3, there appears to be no deviation from Keplerian motion in the vicinity of the bridge. The lack of deviation also means it is likely not an outflow from the inner disk as well.


\section{Conclusion}

In this paper we presented the imagea of the dust continuum at 0.87 mm, the HCO\textsuperscript{+} J=4--3 emission, and the CO J=3--2 emission of the disks around PDS 70. We use these along with the SED of the system and $H$-band polarized data to produce a MCRT model of the system. The main results of our modeling and observation are as follows.

\begin{itemize}

\item There is an azimuthal gap in the HCO\textsuperscript{+} J=4--3 emission which has not been observed in other gas disks. Several possible causes include gravitational depletion by an embedded body, a large stellar flare, chemical changes due to the heating of an embedded body, and shadowing of the gas by an inner disk component. Further study is needed however to determine the most likely candidate.

\item The sub-au region of the disk is gas poor despite evidence for a dusty inner disk component. This suggests that there is likely no gas accretion onto the star.

\item There appears to be an arm in the HCO\textsuperscript{+} J=4--3 emission. The arm could be caused by a planet either very close to the star ($\sim$5 au), or in the radial gap seen in the HCO\textsuperscript{+} J=4--3 and the ring seen the continuum image.

\item The CO J=3--2 data exhibits a bridge along it's minor axis. The bars linear shape and the lack of deviation from Keplerian motion in the  moment 1 data in that region indicate that it lies within the plane of the disk, but has not been smeared out due to Keplerian rotation. The optical thickness of CO gas, and the close proximity of the bridge to the two dim regions in the CO J=3--2 emission, suggest the bridge could be shadowing the gas disk.

\item MCRT modeling, 0.87 mm continuum, and $H$-band polarized data suggest that PDS 70 has a two component dust disk. The inner disk component extends from the sublimation radius out to $\sim$15 au and is followed by a radial gap to $\sim$60 au. The outer disk component then extends from $\sim$60 au to $\sim$110 au. 

\item We do not observe a difference between the small and large dust grain locations in the outer disk as proposed in previous work. Our 0.87 mm continuum image, which probes large grains on the order of 1 mm in size, exhibits a similar radial extent to that of the $H$-band polarized imagery, which probes small grains. The difference between our findings and that of \citet{hashimoto15} is likely due to the low angular resolution of SMA data when compared to ALMA data.

\item The CO J=3--2 moment 1 measurements suggests a dynamical mass of PDS 70 of 0.6 - 0.8 M$\astrosun$ similar to what has been seen in \citet{hashimoto15} and is consistent with the extended kinematics around a K5 star. 
   
\end{itemize}

In this paper we have discussed various features which appear in the gas disk of PDS 70 but are not present in the dust disk. Although we cannot determine a definitive cause for these features with existing data, their existence does provide new incentives to observe this disk with other high angular resolution, high contrast instruments such as VLTI MATISSE or the NIRISS aperture masking interferometry mode of the James Webb Space Telescope. These instruments would allow us to probe the inner regions of the disk where we suspect large planets could be causing many of these features.

We also discussed the overall structure of the dust disk found from continuum data and MCRT modeling. The disk seems to have no strong segregation between the dust grain sizes in the outer disk, but the inner disk does appear to be devoid of large grains. This is contrary to what is suggested in previous work, but is supported by the large silicate peaks in the SED. Further investigation of the inner disk is necessary to confirm its radial extent and optical depth, and may be accomplished through use of the instruments mentioned above.

\acknowledgements

This work is supported supported by the NASA XRP grants NNX17AF88G and NNX16AJ75G. MC thanks the support from the Centro de Astrof\'isica de Valpara\'iso. S.K. acknowledges support from an STFC Rutherford Fellowship (ST/J004030/1) and ERC Starting Grant (Grant Agreement No.\ 639889). This work is supported by the Astrobiology Center Program of National Institutes of Natural Sciences (NINS) (Grant Number: AB281013) and by MEXT KAKENHI No. 17K05399 (E.A.). Y.H. is currently supported by Jet Propulsion Laboratory, California Institute of Technology, under a contract from NASA.

