Revisiting the epoch of the earliest Chinese star catalog titled “Shi Shi Xing Jing”

Abstract

1 Introduction

The observation of stellar positions in Chinese history traces back to 104 BCE, when a famous historian Sima Qian organized a group of astronomers including Luoxia Hong to measure the positions of major reference stars in order to make a calendar system of the Han dynasty’s own. A few decades later, around 77 BCE, another group of astronomers led by Xianyu Wangren made observations to verify the reliability of the calendar system. In 52 BCE, a minister called Geng Shouchang made new observations. These facts were summarized by the philosopher Yang Xiong (53 BCE–18 CE) in his book titled “Fa Yan,” and that Luoxia Hong constructed armillary spheres,¹ Xianyu Wangren made observations with them, and the minister Geng Shouchang made a celestial globe with the observed data. Since these astronomers lived in the Former Han period, it is thought that meridian observations were first performed then (Sun & Kistemaker 1997).

The astronomers belonging to the so-called Shi Shi school attributed their achievements to the legendary figure Shi Shi. Hence, the information on stellar positions was handed down in the name of Shi Shi. Later they were compiled by Chen Zhuo, an astronomer of the Three Kingdoms period (220–280 CE) in Chinese history, into a book titled “Shi Shi Xing Jing” (“Shi Shen’s Astral Canon,” hereafter SSXJ), with supplementary information from the Gan De and Wuxian schools. ² Chen Zhuo compiled 283 asterisms comprising 1464 stars in the Chinese sky.

The book of SSXJ itself was lost, but the coordinates of the reference stars for 118 asterisms including the 28 lunar lodges are preserved in an early Tang astrological compendium called “Kai Yuan Zhan Jing” (ca. 720 CE, hereafter KYZJ) or the “Treatise on Astrology of the Kai Yuan Era.” The coordinate values of the SSXJ stars are partially preserved in other material such as “San Jia Bu Zan” (“Records and Praise for the Three Schools’ Astral Canon”)³ and “Cheonji-Seosangji” (“Treatise on Auspicious Phenomena in Heaven and on Earth”).⁴ The configurations of the asterisms are depicted in the Dunhuang star chart, which dates back to 649–684 CE (Bonnet-Bidaud et al. 2009) and was discovered in 1907 at the Silk Road town of Dunhuang; it is now housed in the British Library. The configurations are also depicted on the Korean ancient star chart called “Cheonsang-Yeolcha-Bunyajido” (Ahn 2011). These data can be used to complement the KYZJ’s coordinate values.

Although SSXJ was named after an astronomer of the fourth century BCE called Shi Shen, its epoch of observation has been a matter of debate. The detailed history of researches on this topic was described in detail by Sun and Kistemaker (1997) and Nakamura (2018). The estimated epoch was first reported by Shinjo to be around 50 BCE, but he later published an article that reported his estimation of the observational epoch to be 360 BCE (Shinjo 1928). His opinions seemed to diverge. Hence, a mathematical verification was performed by Ueda (1929), who concluded that there exist two groups of stars in the determinative stars⁵ of 28 lunar lodges: three of them were observed around 300 BCE, while 18 of them were observed around 200 CE. More recently, Pan (2009) analyzed the whole SSXJ catalog, and reached a similar conclusion that the SSXJ stars are comprised of two groups of stars, observed around 450 BCE and 160 CE, respectively.

Yabuuchi (1937) and Qian (1983) undertook independent textual studies on Chinese classics to conclude that the SSXJ texts were composed during the Former Han dynasty. In particular, Yabuuchi (1937) claimed the epoch of observation to be approximately 70 BCE. Maeyama (1977) pointed out the existence of a misalignment error during the installation of meridian instruments, and performed analytic studies to obtain an epoch of 70 ± 30 BCE.

However, there have been only a couple of research works that gave the statistical uncertainty in the epoch of observation. Recently, Sun and Tian (1993), Sun (1994), and Sun and Kistemaker (1997) also considered the possiblity of meridian observations and a misalignment error. They developed a so-called Fourier analysis method. They analyzed all 120 stars, and concluded that all the stars in the SSXJ catalog were observed around 78 ± 20 BCE. They regarded the standard deviation of the mean for a finite sample of measurements as an estimator of uncertainty for the total population. This quantity will be fine when the number of samples is large, but it can be a rough estimator when the number is small. They did not perform a statistical estimation of the uncertainties by using methods such as bootstrap resamplings. On the other hand, Nakamura (2018) determined the epoch of observation by applying least-squares fitting and bootstrap resamplings to only the 28 determinative stars rather than the 120 stars to obtain the epoch of 54 ± 12 BCE. Moreover, he did not take into account the fact that the polar angles could have misalignment errors caused by measurement with a meridian instrument.

Thus, in this paper we will determine the epoch of observation for all the reference stars listed in SSXJ by using both least-squares fitting and the Fourier method, including instrumental misalignment. In order to estimate the uncertainty, we perform bootstrap samplings. In section 2 we introduce the data and the analysis methods. The results are given in section 3. In section 4 we discuss the statistical implications of the results by performing hypothesis tests.

2 Analysis of the data

2.1 Data

The positions of SSXJ stars are preserved in KYZJ. The original books of KYZJ had disappeared after the Song dynasty, but in 1616 CE a manuscript copy was recovered in the belly of the Buddha’s statue during its restoration (Sun & Kistemaker 1997). This copy was later inserted in Siku Quanshu, whose first workable drafts were completed in 1781 CE. We reference volumes 65–68 of KYZJ in the Wenyuan Library version, which contains the name, the number of constituent stars, a rough description of its location relative to other neighboring asterisms, the lodge angle, the polar angle, and the ecliptic angle for each asterism. For example, the first entry in volume 65 of KYZJ reads:

Shi Shi told that the asterism called Sheti [1] consisting of six stars embraces Dajiao [2] (α Boo). [Its lodge angle is 8 du (du is the Chinese degree) from the reference star of the lunar lodge Jiao [93] (α Vir); its polar angle is 59 and a half du from the north celestial pole; it is at 32.75 du inside of the ecliptic.]

First of all, the lodge angle is the eastward angular difference of a celestial object in right ascension from the reference star of its nearest-to-the-west lunar lodge. Hence, the reference star of Sheti [1] has a right ascension of 8 du plus the right ascension of α Vir. Second, the polar angle is simply the anglular distance from the north celestial pole to an object. We should not be confused that the polar anlges were measured from Polaris. Last of all, the so-called ecliptic angle, called huang dao nei wai du, is the angular distance from the object to the ecliptic along the hour circle that penetrates the object. This is not the ecliptic latitude of the object, which was also pointed out by Yabuuchi (1937) and demonstrated with astronomical calculations by Sun and Kistemaker (1995) and Sun (2000). These ecliptic angles are useful to crosscheck the polar angles.

The identification of stars is the first and crucial step in this analysis. In this paper we basically adopt the identifications made by Sun and Kistemaker (1997). However, we check their identifications one by one. We calculate the equatorial coordinates of the stars at the epoch of 104 BCE by using the PC planetarium software Stellarium. We find a number of exceptional cases, in which other identifications seem to be better. The details are shown in the Appendix. We present those exceptional cases in table 1, where a double dagger at the end of a name represents an exceptional case. A dagger at the end of a name means that the star is discarded during the analysis as an outlier due to its large deviation from the fitting curves.

The lodge angles are given in the seventh column of table 1, along with the reference stars denoted by the numbers of lunar lodges in the sixth column. Here, a number n in the sixth column means the determinative star of the nth lunar lodge. The lodge angles and polar angles are described in Chinese degrees or du, which can be converted into Babylonian degrees by multiplying by 360°/365.25 du. We will use the Babylonian degrees in the main text. We also note that the declination of a reference star, δ_c, is obtained from the polar angle, P, using the relationship δ_c = 90° − P. The right ascensions and the declinations are given in Babylonian degrees. The residuals in declination, Δ_i, are also given in Babylonian degrees for the epoch of 104 BCE.

2.2 Analysis methods

One can find the epoch of observation by finding the equatorial plane of a certain time that fits best to a set of positional data considering the precession of equinoxes. We make use of the fact that the polar angles observed at an epoch t_r deviate from the polar angles at a certain epoch t_s, which is not too far from t_r, form a sinusoidal function of right ascension. The polar angles are assumed to have been measured with meridian instruments when the astronomical objects culminate on the meridian line. Meridian observations have systemic errors because the axis of the meridian instrument may not be aligned with the axis of the celestial poles. As such, we can analyze the observed coordinates of the stars as follows. A more comprehensive description of the methodology is given by Sun and Kistemaker (1997).

Let us denote the deviation of an observed declination |$\delta _i^{\rm obs}$| from the calculated value |$\delta _i^{\rm cal}$| for the ith reference star at a supposed epoch, t_s, by Δ_i. In this paper we will analyze the polar angles of N reference stars, so i = 1, 2, …, N. The given polar angle of the ith reference star, P_i, is converted into the declination |$\delta _i^{\rm obs}$| by |$\delta _i^{\rm obs}=90^{\circ }-P_i$|⁠. We obtain the calculated declination of the star, |$\delta _i^{\rm cal}$|⁠, at the supposed epoch, t_s, by converting the J2000.0 coordinates to the coordinates at the supposed epoch by considering the precession of equinoxes with the method described by Meeus (1998). For the J2000.0 coordinates, we adopt the ICRS J2000.0 coordinate obtained from the SIMBAD.⁶ Although they are small, we also take the proper motions into account. Then, we define the deviation |$\Delta _i\equiv \delta _i^{\rm cal}-\delta _i^{\rm obs}$| for the ith reference star.

These deviations can be fitted with a sinusoidal function f(α) = c + acos α, as suggested by Sun and Kistemaker (1997). Here, α is the right ascension of the reference star at the supposed time of observation t_s. We see that the amplitude a = n × (t_s − t_r), where the rate of change of the declination due to the precession of equinoxes n = 0.°5567 per century, and (t_s − t_r) is the time difference between the epoch of observation t_r and the supposed time t_s. Thus, at the epoch of observation, the fitting parameter a = 0. In other words, in order to find the epoch of observation, we have only to find t_s when a = 0. The coefficient c represents the misalignment angle of the instrument’s axis of rotation with respect to the direction of the north celestial pole in meridian observations.

Thereby, we know that [Δ_i] ≃ f(α_i) = c + acos α_i, where α_i is the right ascension of the ith reference star at the trial epoch t_s with the precessional effects corrected. Here, the coefficients can be obtained by |$a={1\over \pi }\int _{0}^{2\pi }f(\alpha )\cos \alpha \, d\alpha$| and |$c={1\over 2\pi }\int _{0}^{2\pi }f(\alpha )d\alpha$|⁠, which are similar to Fourier analysis. In practice, we can obtain the amplitude a and the misalignment c by using the numerical integration of Δ_i. That is, |$a\simeq {1\over \pi }\sum _{i=1}^{N}[\Delta _i]\cos \alpha _i \, \delta \alpha _i$| and |$c\simeq {1\over 2\pi }\sum _{i=1}^{N}[\Delta _i]\delta \alpha _i$|⁠, where |$\delta \alpha_{i}$| is the interval in right ascension between neighboring reference stars.

However, this numerical scheme is not useful for either highly sparse or highly uneven data sets. Thus, we devise a method to fit the data points with the same function f(α) = c + acos α by minimizing the weighted sum of squared residuals, |$\Phi ^2(c,a)\equiv \sum _{i=1}^{N}w_i r_i^2$|⁠. Here, r_i ≡ Δ_i(α_i) − f(α_i) is the residual for the ith data point, and |$w$|_i is the weight given to that point. The best fit can be found by selecting parameters (c, a) that minimize this quantity Φ²(c, a). Thus, this method is a special case of generalized least squares called weighted least squares (Press et al. 1988; Bevington & Robinson 2003). This can be applied to samples of independent measurements, which is guaranteed by our assumption that all polar angles must have been measured independently with the same meridian instrument on the same occasion.

The weight |$w$|_i for each point can be assumed to be inversely proportional to the variance of each data point. That is to say, |$w_i=1/\sigma _i^2$|⁠. Here, σ_i can be approximated by the measurement error s_i for the ith star. However, no errors are given in the catalog. Thus, we assume that the measurement errors are constant or s_i = s, because they were measured with the same meridian instrument and method on the same occasion. Then, s is the sample standard error estimated from the entire measurement data for the N reference stars. If the number of stars is sufficiently large, σ ≃ s. Therefore, we have only to calculate a sort of merit function |$\phi (c,a)\equiv \Phi ^2(c,a)/w_i=\sum _{i=1}^{N}r_i^2$| to find the parameters (c, a) that minimize φ(c, a) for given data points. Then, we can estimate the measurement error, s, by calculating |$s^2=\sum _{i=1}^{N}r_i^2/(N-1)$| for the best-fitting parameters. If we approximate |$\sigma _i^2\simeq s^2$| or |$w_i=1/\sigma _i^2\simeq 1/s^2$|⁠, then Φ²∝χ². For this reason, this least squares method can also be called a χ² method.

There are a number of extreme outliers whose residual is larger than 2 σ error bars. Hence, we remove them and make another best fit by minimizing their χ² value with the procedure explained above. After a couple of iterations we get sufficiently accurate results. These are the basic procedures for analyzing the coordinates of the SSXJ reference stars.

Note that the uncertainty estimations in the above-mentioned two methods are just approximations. Hence, we make a more robust and sophisticated estimation by adopting a bootstrap method to calculate the confidence interval for the epoch of observation t_r. The number of data points in our case is not large enough to confidently assume that the parent distribution of the residuals is Gaussian. In this case, the bootstrap method is useful for estimating standard errors and obtaining confidence intervals for t_r. We apply Efron’s approach or residual resamplings only to the Δ_i values (Efron & Tibshirani 1993).

The bootstrap samplings are carried out as follows:

1. Make a best fit to the original data with either the Fourier method or the χ² method to find the epoch t_s = t_r.

2. Obtain the residuals r_i = Δ_i − f(α_i) when t_s = t_r.

3. A set of resampled residuals |$\lbrace r_i^*\rbrace$|⁠, i = 1, ⋅⋅⋅, N, of the same size is drawn with replacement from {r_i}, i = 1, ⋅⋅⋅, N.

4. The resampled residuals |$\lbrace r_i^*\rbrace$| are added to the predicted values of the original fit. Here, the predicted values means the predicted declination |$(\delta _i^{\rm obs})^{\prime }=\delta _i^{\rm cal}-f(\alpha _i)$| or |$(\delta _i^{\rm obs})^{\prime }=\delta _i^{\rm cal}-c$| at t_s = t_r, which will yield the best-fitting parameters (c, a) at t_s = t_r. Hence, the resampled declinations are |$\delta _i^{\rm obs*}=(\delta _i^{\rm obs})^{\prime }+r_j$|⁠, where the subscript j is an integer number j = R × N + 1 for a uniform random number 0 < R < 1.

5. Then, make a best fit to the resampled data set |$\lbrace \delta _i^{\rm obs*}\rbrace$| to obtain |$t_{\rm r}^*$|⁠.

6. Repeat steps 3–5 a sufficiently large number of times.

We finally obtain statistics such as mean, median, percentiles, standard error, and confidence intervals for |$\lbrace t_{\rm r}^*\rbrace$|⁠. In this paper we make 40000 bootstrap resamples. The 95% (90%) confidence intervals for t_r can be obtained from the bootstrap results. For example, the 95% (90%) confidence interval [t_A, t_B] can be found by just reading off from the sorted epochs calculated for the resampled data sets. That is to say, t_A is the t_r value for the 2.5 (5) percentile and t_B is the t_r value for the 97.5 (95) percentile of the resultant t_r values. We note that the distributions of |$t_{\rm r}^*$| values in this paper are usually normal. Note that the 1 σ, 2 σ, and 3 σ confidence levels correspond to the |$68.2689\%$|⁠, |$95.4500\%$|⁠, and |$99.7300\%$| confidence levels, respectively.

3 Results

3.1 Entire SSXJ stars

Assuming that both the polar angles and the lodge angles of all the SSXJ stars were measured on the same occasion, we apply the Fourier fit method to all 116 SSXJ stars given by Sun and Kistemaker (1997). As shown in table 2, we obtain the epoch of observation t_r = 〈t_r〉 ± Δt = −79 ± 21, the 1 σ measurement error s = 1.°24, and the misalignment angle c = −0.°82. The uncertainty in the epoch of observation is obtained by generating 40000 bootstrap resamples to obtain the normal distribution of the epoch with |$\langle t_{\rm r}\rangle \pm s/\sqrt{N}=-78.8\pm 33.8\:$|yr (to 1 σ), which is equivalent to the bootstrap 95% confidence interval t_r = [−145, −11] and the 90% confidence interval t_r = [−134, −23]. We show these results in the upper part of table 2.