# Spectroscopic Pulsational Frequency Identification and Mode Determination of
$\gamma$ Doradus Star HD 12901 ††thanks: This paper includes data taken at the
Mount John University Observatory of the University of Canterbury (New
Zealand), the McDonald Observatory of the University of Texas at Austin
(Texas, USA), and the European Southern Observatory at La Silla (Chile).

E. Brunsden1, K. R. Pollard1, P. L. Cottrell1, D. J. Wright2, P. De Cat3
1Department of Physics and Astronomy, University of Canterbury, Private Bag
4800, Christchurch, New Zealand
2Department of Astrophysics, University of New South Wales, Sydney, NSW 2052
Australia
3Royal Observatory of Belgium, Ringlaan 3, 1180 Brussel, Belgium E-mail:
<EMAIL_ADDRESS>
###### Abstract

Using multi-site spectroscopic data collected from three sites, the
frequencies and pulsational modes of the $\gamma$ Doradus star HD 12901 were
identified. A total of six frequencies in the range 1-2 d-1 were observed,
their identifications supported by multiple line-profile measurement
techniques and previously-published photometry. Five frequencies were of
sufficient signal-to-noise for mode identification and all five displayed
similar three-bump standard deviation profiles which were fitted well with
($l$,$m$)=(1,1) modes. These fits had reduced $\chi^{2}$ values of less than
18. We propose that this star is an excellent candidate to test models of non-
radially pulsating $\gamma$ Doradus stars as a result of the presence of
multiple (1,1) modes.

###### keywords:

line: profiles, techniques: spectroscopic, HD12901, stars: variables: general,
stars: oscillations

††pagerange: Spectroscopic Pulsational Frequency Identification and Mode
Determination of $\gamma$ Doradus Star HD 12901 ††thanks: This paper includes
data taken at the Mount John University Observatory of the University of
Canterbury (New Zealand), the McDonald Observatory of the University of Texas
at Austin (Texas, USA), and the European Southern Observatory at La Silla
(Chile).–8††pubyear: 2012

##  1 Introduction

The identification of the geometry of a $\gamma$ Doradus type pulsation
remains one of the more difficult spectroscopic fields. However, it is a
powerful way to characterise the interior structure of a star and improve
stellar pulsation models. The pulsations of $\gamma$ Doradus stars are
particularly of interest as they propagate through the deep layers of a star.
The flux-blocking mechanism at the base of the surface convective zone (Guzik
et al., 2000) is held accountable for the origin of the pulsations which are
g-modes, where the restoring force is gravity.

Briefly $\gamma$ Doradus stars are slightly evolved late-A to late-F stars at
the cooler end of the classical instability strip that display high-order
g-mode pulsations with frequencies on the order of 1 cycle-per-day. For a full
definition see Kaye et al. (1999) and recent reviews of the field are in Kaye
(2007) and Pollard (2009). Currently there are less than 100 bright bona fide
$\gamma$ Doradus stars known (see list in Henry et al. 2011) and a handful of
$\gamma$ Doradus/$\delta$ Scuti hybrid stars (Henry & Fekel, 2005;
Uytterhoeven et al., 2008), with a further 100 $\gamma$ Doradus and 171 hybrid
stars thus far reported by the Kepler mission (Grigahcène et al., 2010;
Uytterhoeven et al., 2011).

Now is truly the age of satellite photometry and no ground-based methods can
compete with the long uninterrupted datasets that satellites such as CoRoT and
Kepler produce. The Fourier frequency spectra of such studies show the high-
levels of precision only obtainable from space (e.g. at least 840 frequencies
found in HD 49434 using the CoRoT satellite; Chapellier et al. 2011). However,
to further our understanding of g-mode pulsations we need more data than
frequencies alone. Several successful techniques using ground-based multi-
colour photometry can be employed to determine the number of nodal lines ($l$)
in a star. The full mode-identification of a star, that is finding the $l$ and
also the number of nodal lines passing through the pole of a star $m$, remains
the sole domain of spectroscopy. It is hoped that modelling of the frequencies
and modes of a star will allow determination of the $n$, number of shells
interior to the stellar surface.

The spectroscopic study of $\gamma$ Doradus stars relies on the collection of
a large amount of high-resolution and preferably multi-site data of sufficient
signal-to-noise for classification of the pulsations. Improvement is gained by
the use of data from multiple sites to reduce daily aliasing patterns in the
Fourier spectra. Such observational campaigns take many months to years in
order to collect a sufficient number of spectra.

Table 1: Observation log of the multi-site data for HD 12901 Observatory | Telescope | Spectrograph | Observation range |  $\Delta$T (d) |  $\\#$ spectra
---|---|---|---|---|---
MJUO | 1.0m | HERCULES | Feb 2009-Sep 2011 | 935 | 478
La Silla | 1.2m | CORALIE | Nov 1998-Nov 2002 | 1166 | 53
McDonald | 2.1m | SES | Sep 2009-Oct 2009 | 8 | 60

When such datasets are compiled for $\gamma$ Doradus stars, the pulsational
frequencies and modes can be examined and compared to those from photometry.
There are still only a handful of $\gamma$ Doradus stars with full mode
identifications and our immediate goal is to classify as many as we can. The
results of these mode identifications can be used to improve pulsational
models (such as those of Grigahcène et al. 2012; Townsend 2003) by providing
information about the amplitudes of excited modes and also to start to
identify patterns within the class.

This paper focuses on the $\gamma$ Doradus star HD 12901 outlining the
observations made and reduction procedure in Section 2. The spectroscopic
frequency analysis for the individual observatories and the combined data
results are described in Section 3 with each method tested. Section 4
describes the reanalysis of white-light and seven-colour photometry taken of
this star. The mode identification of the five identified frequencies follows
in Section 5. Finally Section 6 discusses the findings and their implications
for future work in this field.

##  2 Observations and data treatment

The major findings of this research are the identification of the frequencies
and modes using spectra obtained at three observing sites (Mt John University
Observatory (MJUO), La Silla Observatory and McDonald Observatory) and are
summarised in Table 1.

Spectra were reduced according to the standard local spectrograph software
then they were normalised, continuum fitted and order-merged by the authors
using a semi-automated MATLAB routine. The full spectra were cross-correlated
for each site using a scaled delta function routine (Wright 2008, Wright et
al. 2007). This produces the line profiles used for pulsational analysis. The
signal-to-noise of a cross-correlated line profile is much higher than a
single spectral line and thus pulsations are easier to extract. The line
profiles for each site were combined (scaled, shifted in velocity space and
continuum fitted) to produce a single dataset of consistent line profiles.
This was done post-cross correlation to maximise the number of available lines
used for each spectrum, as different spectrographs operate over different
wavelength ranges.

The line profiles were analysed in FAMIAS, a pulsational frequency mode-
identification toolbox (Zima, 2008), using techniques applied from Zima et al.
(2006). The Fourier spectra for the Pixel-by-Pixel (PbP) technique and that of
the 0th-3rd moments (Balona, 1986; Briquet & Aerts, 2003) were used to
identify frequencies. These frequencies were then analysed to determine mode
identifications using the Fourier Parameter Fit method (Zima, 2009).

Additionally, the software package SIGSPEC (Reegen, 2007) was used to compare
the frequency selection method of FAMIAS. SIGSPEC performs a Fourier analysis
of a two-dimensional dataset and selects frequencies based on their spectral
significance. The frequencies obtained with SIGSPEC can be regarded as more
precise than that of the highest-peak direct selection as, after producing the
Fourier frequency spectrum, the peak with the highest spectral significance is
selected. Spectral significance includes the analysis of the false-alarm
probability to remove frequency peaks caused by irregular data sampling or
noise in the data.

Both FAMIAS and SIGSPEC were used to re-analyse photometric data, originally
published in Handler (1999), Eyer & Aerts (2000) and Aerts et al. (2004), as
further insight into the frequencies of the pulsations. To identify the $l$
values of the modes, the frequency amplitude ratio and phase differences
method, based on pre-computed grids of models, were used (Balona & Stobie,
1979; Watson, 1988; Cugier et al., 1994; Daszyńska-Daszkiewicz et al., 2002).
This was done using the photometric analysis toolbox in FAMIAS (Zima, 2008).

##  3 Spectroscopic Frequency Analysis

In total $591$ spectra from the three sites were of sufficient quality to be
analysed. The data spanned a total of 4667 days, just over 12 years. Each
dataset was analysed independently for frequencies and they were then
combined. By doing this analysis we can comment on the limitations of single-
site data and the extent to which multi-site data reduces aliasing.

###  3.1 Frequency Identification For Each Dataset

####  3.1.1 MJUO

This is the largest single-site dataset and was analysed most extensively.
First the cross-correlated line profiles were analysed in FAMIAS. The
frequency list is in Table 2. The PbP technique found eight possible frequency
peaks including a one-cycle-per-day frequency ($f_{m5}$), likely an alias. The
zeroth moment was too noisy to extract any frequencies. The first and third
moment Fourier peaks each only showed one frequency clearly, $1.186$ d-1 and
$1.184$ d-1 respectively, which is the same frequency given a conservative
uncertainty estimate of $\pm 0.001$ (uncertainties are dealt with more
formally in Section 3.3). The second moment yielded $0.144$ d-1 as the only
viable frequency. Though many more frequencies may be evident in the data,
further extraction of frequencies above the noise was difficult with any
measure of certainty.

The same moments were then tested in SIGSPEC for further analysis. SIGSPEC
found 20 frequencies above a spectral significance level of 5. It is likely
that noise in the data is causing many misidentifications rather than the
possibility that these could all be real, given the noise level in the Fourier
spectra. The frequencies with spectral significance greater than 15 are
reported for moments one to three in Table 2. From the table it is clear that
$f_{m2}$ to $f_{m8}$ (excluding $f_{m5}$) are viable frequencies found using
multiple methods. The addition of further datasets from other observation
sites should improve the signal-to-noise and reduce any aliasing, particularly
1-day aliasing that occurs in the data.

Table 2: Frequencies of HD 12901 found in the MJUO data. Pixel-by-Pixel (PbP) frequencies found using FAMIAS and moments and significances (Sig) found using SIGSPEC are listed. The amplitudes of the PbP frequencies are scaled to the first identified amplitude. Note some frequencies have been identified as $1-f$ aliases of the PbP frequency. The frequency $f_{m5}$ shows identifications near 1 d-1, likely from the data sampling. | PbP | Moments
---|---|---
| Freq. | Amp. | $1^{st}$ | Sig | $2^{nd}$ | Sig | $3^{rd}$ | Sig
| (d-1) | (rel. to $f_{m1}$) | (d-1) |  | (d-1) |  | (d-1) |
$f_{m1}$ | 1.271 | 1 |  |  |  |  |  |
$f_{m2}$ | 1.186 | 1.8 | 1.183 | 24 |  |  | 1.183 | 23
$f_{m3}$ | 1.681 | 1.5 | 0.678 | 26 |  |  | 0.678 | 28
$f_{m4}$ | 1.396 | 1.7 | 0.403 | 20 | 1.396 | 20 |  |
$f_{m5}$ | 1.001 | 1.3 |  |  | 1.001 | 26 |  |
$f_{m6}$ | 1.560 | 1.6 | 0.460 | 16 | 1.560 | 17 | 1.560 | 19
|  |  |  |  |  |  | 0.541 | 18
$f_{m7}$ | 1.216 | 1.6 | 2.214 | 18 | 0.213 | 23 | 2.214 | 20
$f_{m8}$ | 0.244 | 1.1 | 0.756 | 21 |  |  |  |
$f_{m9}$ |  |  |  |  | 0.104 | 16 | 0.899 | 19
$f_{m10}$ |  |  |  |  |  |  | 2.352 | 20
$f_{m11}$ |  |  |  |  |  |  | 0.836 | 17

####  3.1.2 La Silla

The spectra taken on the CORALIE spectrograph have previously been analysed
and published by Aerts et al. (2004) and De Cat et al. (2006). Frequencies of
$1.04\pm 0.28$ d-1 and $1.30\pm 0.30$ d-1 were reported.

Although the small size of the dataset made it difficult to analyse, frequency
peaks were identified at $0.005$ d-1 and $1.26$ d-1 in the PbP measure,
$0.997$ d-1 and $3.01$ d-1 in the first moment and $0.997$ d-1 and $1.999$ d-1
in the third moment. Peaks that occur very close to integer values are less
reliable as they are likely 1-day alias patterns from the observation times
artificially amplified by data sampling effects in the window function. No
identifiable peaks occurred in the zeroth and second moments. Only the $1.26$
d-1 can be identified with any confidence. It is notable that a peak was
apparent in all the above Fourier spectra at $0.23$ d-1.

SIGSPEC was better able to distinguish the 1-day aliases and identified
frequencies of $0.26$ d-1 in the zeroth moment and $0.23$ d-1 in the first and
third moments. These had significances of $3.3$, $3.5$ and $2.7$ respectively,
which are usually regarded as too small to be significant.

####  3.1.3 McDonald

The McDonald data set comprised of 60 observations taken in September and
October 2009. The data were considered to be particularly useful as having
time intervals that overlap some of the MJUO observations, providing some
independent confirmation of the line profile variation. As a stand-alone data
set the shorter wavelength range of the spectrograph, and the larger regions
of telluric lines, reduced the number of stellar lines that could be cross-
correlated from several thousand to around 100 lines. This has a significant
impact on the signal-to-noise of the final line profile. This, and the low
number of observations, meant frequencies detected in this dataset alone were
unconvincing. Despite this, it is noted that frequencies near $1.38$ d-1 and
$1.24$ d-1 are present in the PbP analysis and a peak near $1.67$ d-1 is
visible in the first moment Fourier spectrum.

Figure 1: Scaled line profiles of the three datasets. McDonald (red), La Silla
(blue) and MJUO (black) observations are shown with the mean line profile
(white). All datasets scale well to show a consistent line profile. The data
are provided as an electronic file (see Supporting Information).

###  3.2 Frequency Identification Of Combined Data For Each Method

Figure 2: Effect of the multi-site datasets on the Pixel-by-Pixel Fourier
spectra, the MJUO data on the left and the combined data on the right. Each
series shows the window function for the data at the top then successive pre-
whitening of frequencies $f_{p1}$ to $f_{p5}$. The Fourier spectra for the
Pixel-by-Pixel analysis of the combined data continues to $f_{p10}$. The
frequencies are listed in Table 3.

The cross-correlated line profiles of the three datasets were combined post
cross-correlation to form a single dataset. The line profiles are displayed in
Figure 1. Combining the line profiles instead of the individual spectra was
the best way to preserve the signal-to-noise in the MJUO and La Silla datasets
which had many more lines than the McDonald spectra. There was a small shift
(6 km s-1) in the systemic velocity between the McDonald and the other two
datasets. The McDonald spectra were shifted to ensure consistent data. The
line profiles were scaled to have the same equivalent widths, a parameter
shown to remain constant by the lack of variation in the zeroth moment in the
MJUO data. The Fourier spectra of the MJUO data alone and the combined
datasets are compared in Figure 2.

The spectral windows for the MJUO data alone and all three data sets are very
similar, but, as expected, the addition of the multi-site data reduces the
amplitude of the secondary peaks of the identified frequency. From Figure 2 it
is clear that the addition of the extra data increases the amplitude of the
Fourier peaks. This is an effect of not only the reinforced signal found in
the other sites spectra, but also an artificial amplification due to the
increased base noise level. The frequency peaks remain in nearly the same
positions and show approximately the same repeating pattern. The one-day
aliasing, although still present, has been reduced to make the true peaks more
evident. In general, the frequency identification becomes clearer with the
addition of the multi-site data, despite it being from small sets compared to
the MJUO data.

####  3.2.1 Pixel-by-Pixel (PbP)

The line profiles showed many periodic variations in each pixel. The ten
frequencies found in the combined dataset are listed in Table 3 and the
Fourier spectra displayed in Figure 2. The first two frequencies found were
the $f_{2}$ and $f_{3}$ from Aerts et al. (2004) from photometry and their
$f_{1}$ shows up as our $f_{p4}$. Further frequencies $f_{p3}$ and $f_{p5}$
have credible three bump variations in their line profiles. It is likely that
$f_{p6}$ is a one-day alias frequency arising from the dominance of one single
site (MJUO contributes about $80\%$ of the observations). The frequencies
$f_{p7}$, $f_{p8}$ and $f_{p10}$ showed clear three bump standard deviation
profiles each with smooth phase changes. This increases the likelihood that
they are real frequencies, although the possibility remains that they are from
residual power from a previous poorly-defined frequency. Finally, $f_{p9}$
showed a clear four bump standard deviation profile at $3.14$ d-1 and $4.14$
d-1 frequencies when tested. As none of the previous frequencies showed a
four-bump pattern in the variation it is more likely this is a real frequency
rather than a misidentification of an earlier frequency.

A check on the independence of the frequencies was done by prewhitening the
spectra by 1-day aliases and 1-$f$ aliases of the first five identified
frequencies. Only $f_{p3}\pm 1$ had clear standard deviation profile
variations resembling pulsation and none of the alternate identifications
altered the identification of following frequencies, indicating that they are
all independent.

A second check of the variability of the line profile with different
frequencies involves phasing the data to the proposed frequency and examining
the structure. To see the changes most clearly this was done using the
residuals of the line profiles after subtraction of the mean. The frequencies
$f_{p1}$ to $f_{p5}$ in Figure 3 the residuals on the left side of each plot
show the ’braided rope’ structure typical of non-radial pulsation, but for the
other frequencies this was less clear. If we require pulsations to produce
regular changes in both red and blue halves of the line profile then we must
discard $f_{p6}$, $f_{p7}$ and $f_{p9}$ as candidate frequencies leaving just
$f_{p8}$ and $f_{p10}$ to investigate further.

An investigation into the phase coverage of the data for each frequency can
also be used as an indicator of the reliability of the frequency. All of the
frequencies $f_{p1}-f_{p10}$ were well covered in phase space except $f_{p6}$,
which was poorly sampled between phases 0.2-0.3 and 0.5-0.8 due to the
dominance of a single site in the data.

Table 3: Frequencies from the Pixel-by-Pixel analysis of the combined dataset. The uncertainty estimate for the frequencies is $\pm 0.0002$d-1 as described in Section 3.3. Frequencies with a strike-though were discarded as described in the text. ID | Freq. | Period | Amp. | Phase | Variation
---|---|---|---|---|---
| (d-1) | (d) | (rel. to $f_{p1}$) |  | Explained
$f_{p1}$ | 1.3959 | 0.7164 | 1 | 0.286 | $13\%$
$f_{p2}$ | 1.1863 | 0.8430 | 0.72 | 0.205 | $21\%$
$f_{p3}$ | 1.6812 | 0.5948 | 0.93 | 0.267 | $30\%$
$f_{p4}$ | 1.2157 | 0.8226 | 0.93 | 0.265 | $39\%$
$f_{p5}$ | 1.5596 | 0.6412 | 0.80 | 0.230 | $45\%$
$f_{p6}$ | 1.0004 | 0.9996 | 0.55 | 0.158 | $48\%$
$f_{p7}$ | 1.2465 | 0.8022 | 0.50 | 0.143 | $50\%$
$f_{p8}$ | 1.2743 | 0.7848 | 0.53 | 0.152 | $53\%$
$f_{p9}$ | 3.1392 | 0.3186 | 0.36 | 0.102 | $54\%$
$f_{p10}$ | 2.5803 | 0.3876 | 0.37 | 0.106 | $55\%$

####  3.2.2 Moment Analysis

The zeroth moment analysis shows few peaks that are similar to those from
other methods. There is a cyclic variation around 1-cycle per day but with the
$1.0$ d-1 peak missing due to the combination of multi-site data. The small
amplitude of variations in the second moment indicates we have only small
periodic temperature variations to account for in the line profiles.

The first moment Fourier spectra showed clear peak frequencies similar to the
PbP method. These are shown in the first section of Table 4. The amount of
variation of the first moment was moderately well explained by the selection
of the highest six peaks in the Fourier spectrum. Beyond $f_{f6}$ the
improvements to the explanation of variation are too small to be conclusive.
When the PbP frequencies were extracted from the first moment data (in
approximately the same order as the highest peak frequencies) a better fit to
the variation was found. With the six frequencies, $f_{p1}-f_{p5}$ and
$f_{p8}$, $76\%$ of the variation was removed. The frequencies found using
SIGSPEC closely matched the first four frequencies found in the highest peak
method then reproduced PbP frequencies $f_{p3}$ and $f_{p7}$ as $f_{f9}$ and
$f_{f10}$.

The Fourier spectra of the second moment also showed a few promising
frequencies once the first peak at $0.0006$ d-1 was removed. The frequencies
$f_{s2}$ to $f_{s4}$ match to, or match to one-day aliases of, frequencies
found in the PbP dataset. It is likely that $f_{s6}$ is a misidentification of
$f_{p2}$. The second section of Table 4 catalogues the frequencies found. Both
the highest peak and PbP frequencies accounted for around the same amount of
variation as they identified nearly all the same frequencies.

Frequency peaks found in the third moment are presented in the final section
of Table 4. The first few frequencies found using this method were generally
similar to those found in the PbP method except the double identification of
$f_{t2}$. This is usually the result of small errors in the identification of
a strong frequency, leaving residual power in the Fourier spectrum. It is
likely that this then impacts the following sequence $f_{t5}$ to $f_{t10}$
which all seem close frequencies to other identifications, and possibly one-
day aliases. The PbP frequencies identified explained more variation,
accounting for a total of $83\%$ of the line profile variation.

Overall the PbP variations have been shown to remove the variation from the
data more accurately, demonstrated by the higher percentage of variation
removed for each moment. This leads us to conclude that the frequencies found
using the highest peaks in the moment methods are less reliable. It is notable
that the SIGSPEC identified frequencies are more reliable than the highest
peak selection.

As an additional test of the robustness of the derived frequencies, a
synthetic dataset, modelled on the time spacings and modes identified in
Section 5, was created to compare with the PbP derived frequencies. The
synthetic line profiles were phased and plotted to compare with the real data
in Figure 3. The synthetic profiles strongly match the observed profiles,
strengthening our confidence in the frequency identification.

Table 4: Resulting frequencies from the First, Second and Third Moments found
using peaks in the Fourier spectra and in SIGSPEC

.

First Moment
---
| Highest Peak | SIGSPEC
ID | Freq. | Var. | Freq. | Sig.
$f_{f1}$ | 1.1835 | $16\%$ | 1.18624 | 26
$f_{f2}$ | 1.6100 | $32\%$ |  |
$f_{f3}$ | 1.4017 | $42\%$ | 1.39591 | 26
$f_{f4}$ | 1.2247 | $52\%$ | 1.2156 | 28
| 1.2154 |  | 2.2273 | 21
$f_{f5}$ | 0.6490 | $57\%$ |  | 22
$f_{f6}$ | 1.5627 | $63\%$ | 1.5596 | 20
$f_{f7}$ | 0.7559 | $66\%$ |  |
$f_{f8}$ | 1.5179 | $69\%$ |  |
$f_{f9}$ |  |  | 1.6812 | 22
$f_{f10}$ |  |  | 0.2436 | 22
$f_{f11}$ |  |  | 2.1362 | 17
$f_{f12}$ |  |  | 0.0002 | 17
$f_{f13}$ |  |  | 2.7789 | 15
$f_{f14}$ |  |  | 0.6736 | 15
Second Moment
| Highest Peak | SIGSPEC
ID | Freq. | Var. | Freq. | Sig.
$f_{s1}$ | 0.0006 | $7\%$ |  |
$f_{s2}$ | 0.2159 | $22\%$ |  |
$f_{s3}$ | 0.3767 | $38\%$ |  |
| 1.3962 |  | 1.3959 | 19
$f_{s4}$ | 1.5596 | $46\%$ | 1.5596 | 20
$f_{s6}$ | 2.1786 | $57\%$ |  |
$f_{s7}$ |  |  | 0.1622 | 25
$f_{s8}$ |  |  | 0.9449 | 23
$f_{s9}$ |  |  | 0.9061 | 19
$f_{s10}$ |  |  | 0.7382 | 19
$f_{s11}$ |  |  | 2.6809 | 17
Third Moment
| Highest Peak | SIGSPEC
ID | Freq. | Var. | Freq. | Sig.
$f_{t1}$ | 1.3987 | $17\%$ | 1.3956 | 31
$f_{t2}$ | 1.1864 | $40\%$ | 1.1862 | 25
| 2.1833 |  |  |
$f_{t3}$ | 1.6810 | $49\%$ | 1.6809 | 23
$f_{t4}$ | 1.3489 | $58\%$ |  |
$f_{t5}$ | 1.1661 | $64\%$ |  |
$f_{t6}$ | 2.5624 | $32\%$ | 0.5597 | 21
$f_{t7}$ | 2.7952 | $41\%$ |  |
$f_{t8}$ | 1.1716 | $49\%$ |  |
$f_{t9}$ | 3.5169 | $58\%$ |  |
$f_{t10}$ | 1.2658 | $64\%$ |  |
$f_{t11}$ |  |  | 1.2156 | 27
$f_{t12}$ |  |  | 2.1923 | 20
$f_{t13}$ |  |  | 0.2437 | 18

Figure 3: Line profile residuals phased on the frequencies identified $f_{p1}$
(top) to $f_{p5}$ (bottom) compared with synthetic profiles. The left panel
shows the observed profiles and the right is a synthetic model of the line
profile using the modes identified in Section 5. All are (1,1) modes.

###  3.3 Frequency Results

The frequencies chosen to proceed with in the analysis were those originally
identified in the PbP results, $f_{p1}$ to $f_{p5}$ and $f_{p8}$ as noted in
Table 3. This choice was made based on the prevalence of these frequencies
through the other techniques used, the decreased sensitivity of the PbP
technique to asymmetric variations and generally observed higher signal-to-
noise of the PbP frequencies.

It is clear from the above sections that, due to the inclusion of multi-site
data, 1-day aliasing was not a big problem when identifying frequencies
present in these data.

Identifying the uncertainty range on the frequencies found is a complex task.
When using SIGSPEC we can use the relation of Kallinger et al. (2008):

$\sigma(f)=\frac{1}{T*\sqrt{sig(a)}},$ (1)

where T is the time base of observations in days and $sig(a)$ is the spectral
significance of the frequency. For the frequencies $f_{p1}$ to $f_{p10}$
identified in the first moment method this gives an uncertainty estimate of
$\pm 0.0002$. Additionally we used the estimation from Montgomery & Odonoghue
(1999) who propose

$\sigma(f)=\sqrt{\frac{6}{N}}\frac{1}{\pi{T}}\frac{\sigma(m)}{A}.$ (2)

Here N is the number of observations. The value $\sigma(f)$ is a one-sigma
uncertainty with an amplitude (A) root-mean-square deviation $\sigma(m)$. This
method was used to derive (an underestimate) of the uncertainties of the
pixel-by-pixel frequencies. The results for the frequencies $f_{p1}$ to
$f_{p5}$ ranged from $\pm 0.00013$ to $\pm 0.00018$. Given the observed
differences between the frequencies identified in the PbP and the moment
methods we view $\pm 0.0002$ to be a good estimate of the uncertainty.

The data were also analysed for possible frequency combinations. It was found
that $f_{p5}(1.559)+f_{p3}(1.681)=3.240$d-1 which is possibly the same as
$f_{p7}$(1.246) and strengthens our earlier removal of this frequency.

It is clear from the above sections that many frequencies extracted from the
data are robust as they appear in multiple methods with clear line profile
variations. The first five frequencies in the PbP analysis appear in almost
every analysis method. The high signal-to-noise in the line profiles for these
frequencies makes them suitable for mode identification. The frequency
$f_{p8}$, although showing evidence for being an independent frequency, showed
a misshapen standard deviation profile, which meant there was not enough
certainty in this frequency identification to proceed with a mode
identification. The final PbP frequency $f_{p10}$ showed no evidence of
appearing in other analysis methods and was rejected for mode identification.

Even with all these frequencies identified, variations in the data remain. The
PbP method only removed $50\%$ of the variation and the moment methods up to
$80\%$. Some of the remaining variation is due to the noise in the data and
possible slight misidentifications of the frequencies, but it is likely that
there remain multiple unidentified frequencies below the detection threshold
of the present data. The presence of these further frequencies is expected for
$\gamma$ Doradus stars as they have dozens and sometimes even hundreds of
frequencies identifiable in photometry. The large number of similar amplitude
frequencies make this star challenging to study. However already we are able
to see more frequencies than previously identified in spectroscopy, which
makes this a promising star for further study.

Table 5: Frequencies found in the HIPPARCOS white-light photometry. The first three frequencies match well with $f_{p4}$, $f_{p5}$ and $f_{p7}$.  Frequency (d-1) | Uncertainty | Significance
---|---|---
1.2158 |  $\pm$0.0003 | 7
1.5525 |  $\pm$0.0003 | 6
2.2449 |  $\pm$0.0003 | 6
2.5665 |  $\pm$0.0004 | 5
Table 6: Frequencies (d-1) found in each of the Geneva photometry filters and significances using SIGSPEC. | B | B1 | B2 | G | U | V | V1
---|---|---|---|---|---|---|---
| f. | sig. | f. | sig. | f. | sig. | f. | sig. | f. | sig. | f. | sig. | f. | sig.
$f_{1}$ | 1.215 | 13 | 1.215 | 13 | 1.215 | 12 | 2.218 | 11 | 1.215 | 13 | 1.215 | 11 | 1.215 | 12
$f_{2}$ | 0.09 | 8 | 0.09 | 8 | 0.09 | 8 | 0.393 | 6 | 0.07 | 5 | 0.09 | 7 | 0.095 | 8
$f_{3}$ | 1.396 | 8 | 1.396 | 8 | 1.396 | 8 | 0.071 | 6 |  |  | 1.396 | 7 | 1.396 | 6
$f_{4}$ | 2.567 | 5 | 1.565 | 5 | 3.56 | 4 | 3.32 | 5 |  |  | 4.32 | 4 | 3.32 | 6
$f_{5}$ |  |  |  |  |  |  | 3.19 | 4 |  |  |  |  | 2.95 | 4

##  4 Photometric Frequency and Mode Identification

Two photometric datasets were analysed for frequencies. The first was a time
series of white light observations from the satellite HIPPARCOS (Perryman &
ESA, 1997). Specifically these were taken in the Hp filter. The data span a
period of 1166 days from November 1989 to February 1993 during which 122
measurements were taken. The second dataset was Geneva photometry taken on the
0.7m Swiss telescope at La Silla over a period of 25 years from 1973 to 1997.
A total of 174 observations for each filter in this time. This dataset has
been extensively analysed in Aerts et al. (2004). We present our re-analysis
to complement the spectroscopic results above.

###  4.1 HIPPARCOS Photometric Frequencies

The Fourier analysis of the HIPPARCOS data measurements have a range of 0.09
Hp magnitudes with an average uncertainty of 0.1 Hp magnitudes. It was
difficult to identify clear frequencies from the Fourier spectra, but the
highest peak was at 1.2701 d-1. Using SIGSPEC we were able to identify the
frequencies and their significances as shown in Table 5.

###  4.2 Geneva Photometric Frequencies

The seven filters of the Geneva photometry provide us with enough information
to extract frequencies and identify the $l$ modes. Each of the seven filters
showed very similar variation and this was reflected in the individual Fourier
spectra. Frequency peaks were observed at 1.218 d-1, 1.396 d-1 and 0.0907 d-1
(or one-day aliases) in most filters. To formalise the frequency
identification we used SIGSPEC to identify the frequencies and their
significances as shown in Table 6. These frequencies match well to those
previously found in the same data as is discussed further in Section 6.

###  4.3 Geneva Photometry Modes

The identification of the $l$ value of the spectroscopic frequencies $f_{p1}$
to $f_{p5}$ was attempted using the seven filters in the Geneva photometric
data. The identification was done using the amplitude ratio method in FAMIAS
with modes from $l=1$ to $l=3$ tested in all filters. All five frequencies
were found to be solely consistent with $l=1$ modes distinguished primarily by
the first amplitude ratio. The results are largely identical to Aerts et al.
(2004), Their Figure 6 shows the unique fit to the $l=1$ mode. This method,
although not as powerful as spectroscopic mode identification, gives us
independent support to the $l=1$ mode fits of the star.

##  5 Mode Identification of Spectroscopic Line Profiles

Figure 4: The $95\%$ confidence limit (solid line) of the zero-point fit to
the $v$sin$i$. The dashed lines show the limits at the confidence level and
the dotted line indicates the minimum and hence best fit value.

The mode identification was performed using the PbP frequencies $f_{p1}$ to
$f_{p5}$, as chosen in Section 3.3 ($f_{p8}$ not being considered due to the
distorted line profile variations). The modes of the individual frequencies
were first identified and then a best fit, including all five frequencies, was
computed. Initially the zero-point profile was fitted to determine the basic
line parameters. The best fit is shown in Figure 5(a). The $v$sin$i$ was found
to converge at $63.9$ km s-1 with a $\chi^{2}$ of 103. Figure 4 shows the 95%
confidence limit on the determination of $v$sin$i$ which gives a range of
$63.9\pm 0.5$ km s-1. The zero-point fit set the initial line and stellar
parameters of the space searched. The details of the parameters are given in
Table 8 and the resulting best fit modes are given in Table 9. The unusually
high values for mass and radius are discussed in Section 5.1.

The best fit models found for each individual frequency are plotted in Figures
5(b) \- 5(f). The results in Table 9 and Figures 6(a) to 6(e) show that the
mode identifications are well determined as there is no ambiguity in choosing
the model with the lowest $\chi^{2}$. The inclination appears to lie in the
region $20\degree-45\degree$, with the simultaneous fit giving a value of
$27\degree$ for a $\chi^{2}$ of 16.59. Formally the 95% confidence limit of
this parameter gives $27\degree$ ${}^{+52}_{-12}$, so it is poorly constrained
by the mode identification. Given the $v$sin$i$ of the star, a value near
$30\degree$ is plausible as described further in Section 5.1. The $v$sin$i$
value determined varies slightly between the fits as it is modified to change
the width of the standard deviation profile. The measured $v$sin$i$ of the
star is best determined from the fit to the zero-point profile given above.

The detection of five independent frequencies and the occurrence of a large
number of (1,1) modes makes this star an excellent candidate for further
asteroseismic analysis. The prevalence of the (1,1) modes may indicate
sequencing of the $n$-values, or the number of interior shells. It also is
possible that the (1,1) modes may be linked to the rapid rotation of the star.
This is discussed in the context of all $\gamma$ Doradus stars in Section 6.

###  5.1 Rotation and Pulsation Parameters

The results of a preliminary mode identification were tested using the
rotation and pulsation parameter tools in FAMIAS. Using mass = 1.5 $M_{\odot}$
and radius = 1.7 $R_{\odot}$ (typical for a $\gamma$ Doradus star), $v$sin$i$
= 64 km s-1 and various inclinations indicates the rotational parameters of
the star. Shown in Table 7 are the results for $i=30\degree$, which indicates
that this star is not approaching critical rotational velocity. This is the
case for all tested values of $i$ in the range $i=10\degree-90\degree$. The
rotational frequency of the star is dependent on the inclination, ranging from
4.28 d-1 at $i=10\degree$ to 0.74d-1 at $i=90\degree$. The rotational
frequency for a $\gamma$ Doradus star is expected to be on the same order as
the pulsation frequency and these values fall within this range.

Table 7: Rotational parameters and critical limits for a pulsating star with mass = 1.5 $M_{\odot}$,radius = 1.7 $R_{\odot}$, $v$sin$i$ = 64 km s-1 observed at an inclination of $i=30\degree$. vrot | 128 km s-1
---|---
Trot | 0.67 d
$f_{rot}$ | 1.49 d-1
vcrit | 410 km s-1
$v$sin$i_{crit}$ | 205 km s-1
$i_{crit}$ | 9$\degree$

The pulsation parameter tool can be used to get an indication of the
horizontal-to-vertical pulsation amplitude parameter ($\kappa$) and whether
the pulsation and rotation frequencies fall within the mode determination
limits of FAMIAS. The $\kappa$ value falls between the values 50-800 for
frequencies $f_{1}-f_{2}$, which is above 1.0 as expected for $\gamma$ Doradus
stars. These values then show the natural pulsational frequencies to lie
between 0.2 d-1 and 0.7 d-1. Taking the ratio $f_{rot}$/$f_{co-rot}$ gives
values in the range 1-5 for the determined frequencies and modes (all $m=1$).
FAMIAS has an operational range of $f_{rot}$/$f_{co-rot}\leq 0.5$, being
designed to deal with p-mode pulsations with much larger vertical pulsation
components. The very high values of $\kappa$ for the g-mode pulsations in
$\gamma$ Doradus stars and the high rotation rate of this star mean we are
operating beyond this limit. To obtain a reasonable mode identification the
mass and radius were extended to non-physical values for a $\gamma$ Doradus
star to provide observable amplitudes of pulsation (see Table 8). However the
mass and radius are solely used for the calculation of the ratio $\kappa$, and
do not affect other aspects of the fit. Making these changes allowed us to
consider pulsations with $f_{rot}$/$f_{co-rot}$ values of 0.3-0.5 for the
multi-frequency fit in Table 9.

Including more physics to more accurately describe the effects of rotation of
g-modes in current pulsation models would ultimately solve such problems.
Townsend (2003) investigated the effect of higher rotation on g-modes with
high $n$ and low $l$ by increasing the Coriolis forces. It was seen that the
identification of prograde modes (defined here as $m>0$, in Townsend 2003 as
$m<0$), have smaller differences in the rotating scenarios than retrograde
modes. Smaller values of $m$ are also less affected. There is also a
discussion of this in Wright et al. (2011) who estimate the true limits of the
$f_{rot}$/$f_{co-rot}$ ratio in FAMIAS for various $m$ for a $\gamma$ Doradus
range frequency and found FAMIAS could identify modes up to $f_{rot}$/$f_{co-
rot}=1$ for a $m=1$ mode. These results indicate that the mode identification
is unlikely to be affected by the rotation of this star. We can be confident
that using non-physical values for mass and radius affects only the observed
amplitude of the pulsation and not the mode, itself.

A convincing mode identification of the frequencies $f_{1}-f_{5}$ was
obtained.

Table 8: Stellar parameters used in the mode identification. Fixed values for Teff, [M/H] and $\log g$ were taken from Bruntt et al. (2008) and $v$sin$i$ values from the ranges in Bruntt et al. (2008); De Cat et al. (2006); Dupret et al. (2005); Aerts et al. (2004); Eyer & Aerts (2000).  Parameter | Fixed | Min | Max | Step
---|---|---|---|---
| Value |  |  |
Radius (solar units) |  | 1 | 10 | 0.1
Mass (solar units) |  | 0.5 | 50 | 0.01
Temperature (K) | 7193 |  |  |
Metallicity [M/H] | 0.13 |  |  |
log $g$ | 4.18 |  |  |
Inclination ($\degree$) |  | 0 | 90 | 1
$v$sin$i$ (km s-1) |  | 60 | 75 | 0.1
Table 9: Results of mode identification of all four frequencies individually after a least-squares fit is applied (lsf), and all four frequencies simultaneously (sim). | Mode | $\chi^{2}$ | Inc. | Amp | Ph.
---|---|---|---|---|---
| ID |  | ($\degree$) | (km s-1) |
$f_{1}$ lsf | (1,1) | 7.65 | 43.15 | 0.82 | 0.32
$f_{2}$ lsf | (1,1) | 7.00 | 35.79 | 0.74 | 0.80
$f_{3}$ lsf | (1,1) | 7.57 | 22.40 | 1.50 | 0.18
$f_{4}$ lsf | (1,1) | 10.43 | 33.11 | 1.50 | 0.72
$f_{5}$ lsf | (1,1) | 5.49 | 33.11 | 0.50 | 0.61
$f_{1}$ sim | (1,1) | 16.59 | 29.67 | 2.06 | 0.31
$f_{2}$ sim | (1,1) |  |  | 1.07 | 0.81
$f_{3}$ sim | (1,1) |  |  | 2.90 | 0.18
$f_{4}$ sim | (1,1) |  |  | 1.07 | 0.71
$f_{5}$ sim | (1,1) |  |  | 2.17 | 0.60

(a) Fit of the zero-point profile ($\chi^{2}=104$).

(b) Mode identification of $f_{p1}=1.3959$ d-1 ($\chi^{2}=7.7$).

(c) Mode identification of $f_{p2}=1.1863$ d-1 ($\chi^{2}=7.0$).

(d) Mode identification of $f_{p3}=1.6812$ d-1 ($\chi^{2}=7.6$).

(e) Mode identification of $f_{p4}=1.2157$ d-1 ($\chi^{2}=10.4$).

(f) Mode identification of $f_{p5}=1.5596$ d-1 ($\chi^{2}=5.5$).

Figure 5: The fit (dashed) of the mode identification to the mean line
profile, variation and phase (solid) of the five identified frequencies. All
have been identified as (1,1) modes with the labelled $\chi^{2}$ for the best
fit. An indication of the maximum uncertainty is given on each plot.

##  6 Discussion

This study allowed a direct comparison between single-site and multi-site
data. From this we can judge the usefulness of large single-site datasets. The
findings show that although the addition of further sites increased the
amplitudes of the frequencies, it also elevated the base noise level.
Additionally it is clear from the window function spectra in Figure 2 that the
1-day aliasing pattern is reduced but not entirely eliminated. The above leads
us to conclude that the addition of multi-site data is useful, but not
required to extract frequencies from sufficiently large datasets. We must also
require that the data from any one site have a sufficient number of
observations to produce a balanced (well sampled in phase space) mean line
profile in order to combine the cross-correlated line profiles with the
highest precision.

The frequencies found in the PbP method were the most reliable and were
consistent with all the other analysis methods. Past papers have analysed the
photometric data from HIPPARCOS and also some multi-coloured photometry. In
Handler (1999), the author finds a frequencies of 2.18 d-1 and 1.2155 d-1
($f_{p4}$). In this paper we found frequencies $f_{p4}$, $f_{p5}$ and (2x)
$f_{p8}$ in the HIPPARCOS photometry.

The multi-colour photometry was previously published by Aerts et al. (2004),
who found frequencies equivalent to $f_{p4}$, $f_{p1}$ and a 1-day alias of
$f_{p3}$. In the same data we find $f_{p4}$, $f_{p1}$ and $f_{p5}$. The mode
identification showed the frequencies $f_{p4}$, $f_{p1}$ and $f_{p3}$ to best
fit $l=1$ modes, the same as found in this analysis. The same result was found
by Dupret et al. (2005) who showed that the three frequencies from Aerts et
al. (2004) are $l=1$ modes when using time-dependent convection models.

Aerts et al. (2004) also describe the spectroscopic dataset taken with CORALIE
that was also used in this analysis. The authors did not find any frequencies
from the spectroscopic data alone (it is, as the authors note, too small a
dataset for spectroscopic analysis). The imposition of the $f_{p4}$, $f_{p1}$
frequencies found in their photometric analysis did provide some harmonic fits
with low amplitudes.

(a) $\chi^{2}$ values for $l=0-3$ for $f_{p1}=1.3959$d-1.
Best fit value is for ($l,m$) = (1,1).

(b) $\chi^{2}$ values for $l=0-3$ for $f_{p2}=1.1863$d-1.
Best fit value is for ($l,m$) = (1,1).

(c) $\chi^{2}$ values for $l=0-4$ for $f_{p3}=1.6812$d-1.
Best fit value is for ($l,m$) = (1,1).

(d) $\chi^{2}$ values for $l=0-3$ for $f_{p4}=1.2157$d-1.
Best fit value is for ($l,m$) = (1,1).

(e) $\chi^{2}$ values for $l=0-3$ for $f_{p5}=1.5596$d-1.
Best fit value is for ($l,m$) = (1,1).

Figure 6: Lowest $\chi^{2}$ values for each possible ($l,m$) combination for
the final identified frequencies. The best fit $\chi^{2}$ is identified in
bold.

A study applying the frequency ratio method to the frequencies found in Aerts
et al. (2004) was done by Moya et al. (2005), who found three models
consistent with $l=1$ modes. The models have a $T$= $6760K$, log$g=3.88-4.12$
and stellar ages around $2-3$ Myr.

Considering all the prior studies, all of the candidate frequencies we confirm
($f_{p1}$ to $f_{p5}$ and $f_{p8}$) are well supported. The mode
identification is also partially confirmed by photometry. There is no
suggestion of different frequencies observed in photometry and spectroscopy as
there are from some $\gamma$ Doradus stars (see Maisonneuve et al. 2011;
Uytterhoeven et al. 2008).

A prevalence of (1,1) modes in this star, and in $\gamma$ Doradus stars in
general, is beginning to emerge. This includes two modes in HD 135825
(Brunsden et al., 2012), two in $\gamma$ Doradus (Balona et al., 1996; Dupret
et al., 2005), one in HD 40745 (Maisonneuve et al., 2011), one in HR8799
(Wright et al., 2011), two in HD 189631 (Davie, 2011) and two in HD 65526
(Greenwood, 2012). This is possibly due to the large surface area covered in
each segment of a (1,1) mode, meaning pulsations have larger amplitudes and
large changes in amplitude across the stellar surface. This indicates an
observational selection effect. Additionally Balona et al. (2011) discuss the
prevalence of $m=\pm 1$ modes in a sample of $\gamma$ Doradus candidates from
Kepler photometry. These authors assume the dominant frequency to be the
rotational frequency and they then show further frequencies to be close to
this value. This constrains the light maxima to once per rotation cycle,
requiring a $m=\pm 1$ mode. The occurrence of five of these such modes in this
star suggests some physical linking between them. The identification of
several (1,1) modes in this star led to an investigation into the period
spacings of the six identified frequencies (Table 3). The asymptotics of
oscillation theory (Tassoul, 1980) predicts a characteristic period spacing
for high-order g-modes of the same low degree ($l$) for sequential values of
$n$. An investigation into the period spacings of the PbP identified
frequencies shows that the spacing between $f_{p1}$-$f_{p2}$ and
$f_{p1}$-$f_{p3}$ to be close ($0.1266$d and $0.1215$d respectively). This
suggests they could be subsequent values of $n$ if we allow our frequencies to
vary by $\pm 0.003$ d-1. The frequencies $f_{p4}$, $f_{p5}$ and $f_{p8}$
however do not fit with this spacing and additionally it may be that the close
spacing of $f_{p2}$ and $f_{p4}$ is inconsistent with current theoretical
models. The identification of the frequencies could be further improved by
photometric studies to confirm this. Ultimately the sequencing of $n$-values
could provide us with direct information about the stellar interior.

For the reasons above we propose this star as an excellent candidate with
which to test asteroseismic models and potentially give us insight into the
pulsational behaviour of $\gamma$ Doradus stars.

##  7 Acknowledgements

This work was supported by the Marsden Fund administered by the Royal Society
of New Zealand.

The authors acknowledge the assistance of staff at Mt John University
Observatory, a research station of the University of Canterbury.

We appreciate the time allocated at other facilities for multi-site campaigns,
particularly McDonald Observatory and La Silla (European Southern
Observatory).

Gratitude must be extended to the numerous observers who make acquisition of
large datasets possible. We thank P. M. Kilmartin at MJUO and all the
observers at La Silla and the HIPPARCOS team for their dedication to acquiring
precise data.

This research has made use of the SIMBAD astronomical database operated at the
CDS in Strasbourg, France.

Mode identification results obtained with the software package FAMIAS
developed in the framework of the FP6 European Coordination action HELAS
(http://www.helas-eu.org/).

We thank our reviewer Gerald Handler for his helpful comments that improved
this paper.

##  8 Supporting Information

Additional Supporting information may be found in the online version of this
article:

Data file. The line profile data including Modified Julian Date, velocity (on
a relative scale) and the intensity at each of the 120 velocity sampling
points for each profile. Line 1 is the axis scale, lines 2-478 are MJUO
observations, lines 479-526 are La Silla observations and lines 527-586 are
McDonald observations.

## References

  * Aerts et al. (2004) Aerts C., Cuypers J., De Cat P., Dupret M. A., De Ridder J., Eyer L., Scuflaire R., Waelkens C., 2004, A&A, 415, 1079
  * Balona (1986) Balona L. A., 1986, MNRAS, 219, 111
  * Balona et al. (1996) Balona L. A., Böhm T., Foing B. H., Ghosh K. K., Janot-Pacheco E., Krisciunas K., Lagrange A.-M., Lawson W. A., James S. D., Baudrand J., Catala C., Dreux M., Felenbok P., Hearnshaw J. B., 1996, MNRAS, 281, 1315
  * Balona et al. (2011) Balona L. A., Guzik J. A., Uytterhoeven K., Smith J. C., Tenenbaum P., Twicken J. D., 2011, MNRAS, 415, 3531
  * Balona & Stobie (1979) Balona L. A., Stobie R. S., 1979, MNRAS, 189, 649
  * Briquet & Aerts (2003) Briquet M., Aerts C., 2003, A&A, 398, 687
  * Brunsden et al. (2012) Brunsden E., Pollard K. R., Cottrell P. L., Wright D. J., De Cat P., Kilmartin P. M., 2012, MNRAS, 422, 3535
  * Bruntt et al. (2008) Bruntt H., De Cat P., Aerts C., 2008, A&A, 478, 487
  * Chapellier et al. (2011) Chapellier E., Rodríguez E., Auvergne M., Uytterhoeven K., Mathias P., Bouabid M.-P., et al., 2011, A&A, 525, A23
  * Cugier et al. (1994) Cugier H., Dziembowski W. A., Pamyatnykh A. A., 1994, A&A, 291, 143
  * Daszyńska-Daszkiewicz et al. (2002) Daszyńska-Daszkiewicz J., Dziembowski W. A., Pamyatnykh A. A., Goupil M.-J., 2002, A&A, 392, 151
  * Davie (2011) Davie M., 2011, Astr480 project, University of Canterbury
  * De Cat et al. (2006) De Cat P., Eyer L., Cuypers J., Aerts C., Vandenbussche B., Uytterhoeven K., Reyniers K., Kolenberg K., Groenewegen M., Raskin G., Maas T., Jankov S., 2006, A&A, 449, 281
  * Dupret et al. (2005) Dupret M.-A., Grigahcène A., Garrido R., De Ridder J., Scuflaire R., Gabriel M., 2005, MNRAS, 360, 1143
  * Eyer & Aerts (2000) Eyer L., Aerts C., 2000, A&A, 361, 201
  * Greenwood (2012) Greenwood A., 2012, Astr391 project, University of Canterbury
  * Grigahcène et al. (2010) Grigahcène A., Antoci V., Balona L., Catanzaro G., Daszyńska-Daszkiewicz J., Guzik J. A., Handler G., Houdek G., Kurtz D. W., Marconi M., Monteiro M. J. P. F. G., et al., 2010, ApJ, 713, L192
  * Grigahcène et al. (2012) Grigahcène A., Dupret M.-A., Sousa S. G., Monteiro M. J. P. F. G., Garrido R., Scuflaire R., Gabriel M., 2012, MNRAS, 422, L43
  * Guzik et al. (2000) Guzik J. A., Kaye A. B., Bradley P. A., Cox A. N., Neuforge C., 2000, ApJ, 542, L57
  * Handler (1999) Handler G., 1999, MNRAS, 309, L19
  * Henry & Fekel (2005) Henry G. W., Fekel F. C., 2005, AJ, 129, 2026
  * Henry et al. (2011) Henry G. W., Fekel F. C., Henry S. M., 2011, AJ, 142, 39
  * Kallinger et al. (2008) Kallinger T., Reegen P., Weiss W. W., 2008, A&A, 481, 571
  * Kaye (2007) Kaye A. B., 2007, Communications in Asteroseismology, 150, 91
  * Kaye et al. (1999) Kaye A. B., Handler G., Krisciunas K., Poretti E., Zerbi F. M., 1999, PASP, 111, 840
  * Maisonneuve et al. (2011) Maisonneuve F., Pollard K. R., Cottrell P. L., Wright D. J., De Cat P., Mantegazza L., Kilmartin P. M., Suárez J. C., Rainer M., Poretti E., 2011, MNRAS, 415, 2977
  * Montgomery & Odonoghue (1999) Montgomery M. H., Odonoghue D., 1999, Delta Scuti Star Newsletter, 13, 28
  * Moya et al. (2005) Moya A., Suárez J. C., Amado P. J., Martin-Ruíz S., Garrido R., 2005, A&A, 432, 189
  * Perryman & ESA (1997) Perryman M. A. C., ESA eds, 1997, The HIPPARCOS and TYCHO catalogues. Astrometric and photometric star catalogues derived from the ESA HIPPARCOS Space Astrometry Mission Vol. 1200 of ESA Special Publication
  * Pollard (2009) Pollard K. R., 2009, in J. A. Guzik & P. A. Bradley ed., American Institute of Physics Conference Series Vol. 1170 of American Institute of Physics Conference Series, A Review Of $\gamma$ Doradus Variables. pp 455–466
  * Reegen (2007) Reegen P., 2007, A&A, 467, 1353
  * Tassoul (1980) Tassoul M., 1980, ApJS, 43, 469
  * Townsend (2003) Townsend R. H. D., 2003, MNRAS, 343, 125
  * Uytterhoeven et al. (2008) Uytterhoeven K., Mathias P., Poretti E., Rainer M., Martín-Ruiz S., Rodríguez E.and Amado P. J., Le Contel D., Jankov S., Niemczura E., Pollard K. R., Brunsden E., Paparó M., Costa V., et al., 2008, A&A, 489, 1213
  * Uytterhoeven et al. (2011) Uytterhoeven K., Moya A., Grigahcène A., Guzik J. A., Gutiérrez-Soto J., Smalley B., Handler G., Balona L. A., Niemczura E., Fox Machado L., Benatti S., Chapellier E., Tkachenko A., et al., 2011, AAP, 534, A125
  * Watson (1988) Watson R. D., 1988, APSS, 140, 255
  * Wright (2008) Wright D. J., 2008, PhD thesis, University of Canterbury
  * Wright et al. (2011) Wright D. J., Chené A.-N., De Cat P., Marois C., Mathias P., Macintosh B., Isaacs J., Lehmann H., Hartmann M., 2011, ApJL, 728, L20
  * Wright et al. (2007) Wright D. J., Pollard K. R., Cottrell P. L., 2007, CoAst, 150, 135
  * Zima (2008) Zima W., 2008, CoAst, 157, 387
  * Zima (2009) Zima W., 2009, A&A, 497, 827
  * Zima et al. (2006) Zima W., Wright D., Bentley J., Cottrell P. L., Heiter U., Mathias P., Poretti E., Lehmann H., Montemayor T. J., Breger M., 2006, A&A, 455, 235


# A Novel Reconfigurable Vector-Processed Interleaving Algorithm for a DVB-
RCS2 Turbo Encoder

Ohad Boxerman1, Moshe Bensimon1, Shlomo Greenberg and Yehuda Ben-Shimol 1 M.
Bensimon, and O. Boxerman Contribute Equally.

###### Abstract

Turbo-Codes (TC) are a family of convolutional codes enabling Forward-Error-
Correction (FEC) while approaching the theoretical limit of channel capacity
predicted by Shannon’s theorem. One of the bottlenecks of a Turbo Encoder (TE)
lies in the non-uniform interleaving stage. Interleaving algorithms require
stalling the input vector bits before the bit rearrangement causing a delay in
the overall process. This paper presents performance enhancement via a
parallel algorithm for the interleaving stage of a Turbo Encoder application
compliant with the DVB-RCS2 standard. The algorithm efficiently implements the
interleaving operation while utilizing attributes of a given DSP. We will
discuss and compare a serial model for the TE, with the presented parallel
processed algorithm. Results showed a speed-up factor of up to 3.4 Total-
Cycles, 4.8 Write and 7.3 Read.

###### Index Terms:

Digital Signal Processing, Digital Video Broadcasting-Return Channel
Satellite, Permutations, Turbo Codes, Vector Processor, Very Large Instruction
Word.

##  I Introduction

Modern age technology and the use of communicative devices dictate the need
for efficient communication protocols, data transfers and eminent data
manipulation. Nowadays, applications, such as video broadcasting and satellite
communications, rely upon fast encoding and calculation tools. Parallelization
of data processing requires, in some cases, innovative thinking and different
memory access routines unlike the common-default ordered memory structures
used. These approaches are the cause of major performance enhancements of many
applications. Modern processors use multiple cores for processing data, but in
many cases, the program that implements the processing algorithm is sequential
and unsuited for the multi-core-vector attributes of the processor, hence not
exploiting the full potential of the processor. For these reasons, there is a
need for generic tools and techniques for converting known sequential
algorithms to a more efficient parallel implementation over vector and multi-
core processors existing in today’s market.

###  I-A Error Correction

Wireless channels impose multiple noise deficiencies that affect the
transmitted information. High data transfer rates require fast and efficient
encoding that impose both rapid data transmission and eminent data integrity.
Data integrity is achieved by data encoding that enables both error detection
and correction at the receiver, that is, correction of the received data
without the need of resending damaged data packets.

TC implements an approach to the mathematical bounds of data rates with
respect to the given noise over given channels, found by C.E. Shannon [1, 2].
TC adds redundant bits to the original message to enable Forward Error
Correction (FEC), detection and correction of errors at the receiving-end,
thus improving data integrity. These techniques are based on the fact that
some communication channels suffer from noises over certain segments of the
frequency spectrum, and adding redundant bits with a known permutation enables
a reliable detection and correction of defected messages with high
probability. Additional examples of FEC algorithms include Reed-Solomon codes,
Low-Density Parity-Check, Repeat-Accumulate codes and Product codes [3].

###  I-B Turbo Encoders

TC, a subdivision of the Parallel Concatenated Convolutional Codes family are
widely used, mainly because of their adjust-ability for the purpose of real-
time implementations (e.g. in satellite communications [4]). TC are in use in
wireless ATMs, Third-Generation systems and in video broadcasting [5].
Discussion on the topic could be found in [6, 7, 8]. Every TE is constructed
by two main operation blocks: an interleaver and a convolution-encoder. While
the encoding is the purpose of a TE, interleaving causes substantial operation
latency and performance skew. Interleaving is the action of creating a
permutation of the input data by a deterministic shuffle. Modulo calculations
are a common tool used in programming implementations in order to keep the
data in line with a finite vector size or memory allocated array. Modulo
calculations are a hardware obstacle which needs to be addressed, especially
when dealing with various and large modulo bases. In the TE interleaving
block, modulo is a repeated operation consuming both power and computation
time. This article presents a vector-oriented implementation of the
permutation stage of a specific TE described by the DVBRCS2 standard. Sec. II
elaborates the specifications of the turbo code implemented in this research.

###  I-C Vector Processor and Parallel Processing

A Vector Processor (VP) differs from a Scalar Processor (SP) by its ability to
process a single instruction on multiple data (SIMD), usually a one-
dimensional array (vector). That 2 is, a single operation is performed
simultaneously on N independent elements (N being the given length of the
vector). The VP can fetch N elements using a single load operation, thus
saving time in both fetching and decoding the data and instructions. In a VP’s
architecture it is essential to redefine operations while indicating whether
the operation is performed on a scalar or a vector of certain size. The
vectors could be of constant/dynamic length. Amongst the applications suitable
that can benefit from the accelerated performance on a VP are multimedia,
signal processing, cryptography etc. Since the interleaving stage of the TE is
a manipulation of long bit-arrays, the choice of a VP is a natural one.
Requiring the entire packet of bits for the interleaving stage causes it to be
a bottleneck in the TE application. Our simulations revealed that 78% of
processing cycles were of this stage, therefore, improving the performance of
the permutation algorithm will result with speedup of the whole application.
Using a VP allows to efficiently manipulate the input data by adequate vector
and parallel operations, thus enhancing performance of the runtime application
both in memory access (R/W) and in total cycles. The algorithmic solution
presented in this paper is general and can be reconfigured by the basic DVB-
RCS2 parameters.

###  I-D Related Work

With growing demand for fast data transfers, in recent years we notice an
increase in hardware accelerators specifically designed for turbo encoding. A
number of hardware accelerators exists: (1) Texas Instruments offers the
TCI6618 Multicore SoC DSP with up to 582Mbps throughput. (2) NXP
Semiconductors developed the MSC8157 SoC DSP which can reach up to 330 Mbps
throughput [9]. Algorithms matching vector-abilities and smart SIMD processing
with given VP can achieve processing speedup of up to 35% [10] while more
recent work of comparing vectorized vs. non-vectorized execution achieved up
to 57.71% performance improvement [11]. This article reveals a novel algorithm
solving the modulo calculation problem by exploiting the advantages of a VP,
thus avoiding the modulo calculation itself. The organization of this paper is
as follows: Sec. II, describes the specific TE implemented set by the DVB-RCS2
standard; Sec. III, describes the original serial permutation algorithm and
implementation. This algorithm serves as a reference for the speedup factor
calculations; Sec. IV, describes the proposed vectorial and parallel
permutation algorithm; Sec. V, describes the simulation-based analysis and
depicts the results graphically; Finally, Sec. VI, concludes the results and
achievements and proposes future research possibilities.

##  II The Turbo Encoder

Fig. 1 depicts the TE implemented by the DVB-RCS2 standard. It has a coding
rate of 1:3 meaning for every input bit, there are 3 output bits. The input is
a 2-bit stream denoted A and B. The output codeword is combined of 3 couplets:
(1) The original A and B couple unchanged. (2) A and B are encoded by the
‘Encoder Core’, a convolutional calculation block. (3) N couplets are delayed
and rearranged in the ‘Permutation’ block, then encoded via the ‘Encoder Core’
block.

The permutation stage of the TE is determined by five parameters,
$P,Q_{0},Q_{1},Q_{2}$ and $Q_{3}$, with ranges defined by the DVB-RCS2
standard, and vector size N, the number of couplets in bits. The specific
parameters used in the implementation are in compliance with the DVB-RCS2
standard and are detailed in Table I[4].

  1. 1.

The original A and B couple unchanged.

  2. 2.

A and B are encoded by the ‘Encoder Core’, a convolutional calculation block.

  3. 3.

$N$ couplets are delayed and rearranged in the ‘Permutation’ block, then
encoded via the ‘Encoder Core’ block.

Figure 1: Turbo Encoder Scheme

The permutation stage of the TE is determined by five parameters,
$P,Q_{0},Q_{1},Q_{2}$ and $Q_{3}$, with ranges defined by the DVB-RCS2
standard, and vector size $N$ (the number of couplets in bits). The specific
parameters used in the implementation are in compliance with the DVB-RCS2
standard and are detailed in Table I[4].

TABLE I: Example of used Turbo Encoder parameter-sets N[bits] | P | $Q_{0}$ | $Q_{1}$ | $Q_{2}$ | $Q_{3}$
---|---|---|---|---|---
56 | 9 | 2 | 2 | 8 | 0
152 | 17 | 9 | 5 | 14 | 1
236 | 23 | 10 | 2 | 11 | 1
384 | 25 | 1 | 2 | 0 | 1
432 | 29 | 1 | 4 | 1 | 1
492 | 31 | 0 | 3 | 1 | 0
520 | 31 | 0 | 1 | 2 | 0
776 | 39 | 7 | 0 | 0 | 0
1056 | 43 | 0 | 0 | 6 | 2
1192 | 49 | 0 | 3 | 5 | 0
2396 | 81 | 1 | 2 | 5 | 2

Initial simulations running a simple existing serial implementation of this TE
resulted in the permutation stage taking up to 78% of the application (cycle-
wise) and showed potential parallelism features such as repetitiveness modulo
calculations which could be reconfigured to enhance performance. Fig. 2
depicts an example of permuted indexes by original indexes (for $N=776$ and
matching parameters from Table I) and shows the constant $4P$ incrimination
leading to the new parallel approach described in Sec. IV. The turbo encoder
permutation stage is carried out in two levels:

Figure 2: Permuted indexes for example vector length of 776. $\Delta=4\cdot P$
is the offset between two consecutive indexes after permutation.

(1) describes a swap of the bit couplets for every odd indexed component

$\mathrm{if}\\{j\ (\mathrm{mod}\
2)=1\\}\Rightarrow(A_{j},B_{j})\rightarrow(B_{j},A_{j}),$ (1)

(2) describes a calculated rearrangement of the indexes based on the chosen
parameters. Further elaboration of the latter can be found in Section IV.

$i\equiv\Pi(j)=(P\cdot j+Q+3)\ (\mathrm{mod}\ N)$ (2)

where

$Q=\begin{cases*}0,&$j\ (\mathrm{mod}\ 4)=0$\\\ 4Q_{1},&$j\ (\mathrm{mod}\
4)=1$\\\ 4Q_{0}\cdot P+4Q_{2},&$j\ (\mathrm{mod}\ 4)=2$\\\ 4Q_{0}\cdot
P+4Q_{3},&$j\ (\mathrm{mod}\ 4)=3$\\\ \end{cases*}$ (3)

##  III Serial Implementation

The serial implementation serves as reference for speedup and relies on Look
Up Tables (LUT). For every input vector, given its length and matching
parameters (taken from Table I), the program (i) creates LUT by deterministic
calculations matching the permutation indexes of each input bit, (ii) loads
each bit and its designated new index and (iii) stores the bit in a new output
vector. This process is ineffective for two main reasons:

  1. 1.

The LUT calculations are determined the DVB-RCS2 standard. The Modulo based
calculations are complex for any processor and is a changing parameter
disabling simple hardware solutions for these dynamic calculations. Given a
finite collection of sets of vector length and parameters, saving all pre-
calculated LUT is an option. However, real-time processors run on finite and
usually small sized program memory space (MS) bounding the MS that can be
reserved for LUT.s

  2. 2.

Saving the pre-calculated LUT in the memory doesn’t solve the following issue.
Once the LUT are calculated, the program loads bit-by-bit with it’s matching
new index from the LUT then the bit is stored to it’s designated index in the
output vector. Loading and storing bit-bybit (and index-by-index) results in
excessive memory access.

The presented vector algorithm solves these problems by processing the data in
max sized streams (loading and storing more than one bit at a time) and by
replacing the modulo based calculations of the LUT with data and memory
manipulations resulting with lower run-time and memory access. The asymptotic
tight bound computational complexity of the serial algorithm is $O(N2)$.

##  IV New Approach: Vector Algorithm

Using a Vector DSP with Very Large Instruction Word (VLIW) architecture
requires a different approach compared to the original straightforward
approach detailed in Sec. III. The main idea is based on manipulating the data
vectorwise instead of bit-by-bit and taking advantage of the modulo attributes
originating from (2) and (3). Fig. 2 indicates modulo induced strides from one
index to the next. There are four different groups of indexes derived from the
four cases in (3). These strides constitute the basis of the developed
vectorwise algorithm described in the following subsections. The algorithm is
based on configurable parameters defined by a standard ensuring a
comprehensive and generic solution, thus enabling further research
possibilities elaborated in Sec. VI.

Basing the algorithm on Load-Execute-Store operations fitted with the CEVA-
XC4500 DSP attributes creates phases. Each phase consists of iterations of
execution operations that are done in one work-cycle of the DSP. The following
subsections elaborate these operations which combined implement the
permutation block of the TE. The presented algorithm was constructed for a
specific VP and a specific TE but could be easily modified to various VPs and
or different TEs. We chose to test and demonstrate the method on the CEVA-
XC4500 and implement the TE described by the DVB-RCS2 standard.

###  IV-A Bit to Word

As in other DSP architectures, the CEVA-XC4500 DSP requires padding of the
original input-data for more efficient data manipulation. Hence, the first
phase in the algorithm is padding every input-data bit with 15 zeros.
Consequently, every original input-data bit is now represented as one Word (16
bits), as shown in Fig. 3.

Figure 3: Bit to Word - Padding every bit with 15 zeros.

The CEVA-XC4500 DSP is capable of storing 512 bits in one store operation.
Following that $16\cdot 32=512$, loading 32 bits of the original input-data
stream will resolve in 512 bits of padded data as the result of one iteration
of this phase. The number of iterations needed for all the original inputdata
stream to be processed $\big{\lceil}\frac{2N}{32}\big{\rceil}$ , where N is
the number of couplets in bits (see Fig. 1). Hence, the output vector of this
phase is of size $16\times 2N[bits]$ or $2N[words]$. Notes:

The CEVA-XC4500 DSP is capable of storing 512 bits in one store operation.
Following that $16\cdot 32=512$, loading 32 bits of the original input-data
stream results with 512 bits of padded data in one iteration of this phase.
Therefore, the number of iterations needed for all the original input-data
stream to be processed is $\big{\lceil}\frac{2N}{32}\big{\rceil}$, where $N$
is the number of couplets in bits (see Fig.1). Hence, the output vector of
this phase is of size $16\times 2N[bits]$ or $2N[words]$.

Notes:

  1. 1.

The permutation stage is identical to every two bits of original input couplet
(A and B) and so we refer to every padded couplet as a Double-Word (DW). In
addition, The architecture of the CEVA-4500 DSP is designed to processes DWs,
therefore the next two phases of permutation are performed on $N$ DWs.

  2. 2.

The Bit-to-Word phase is added as an input rearrangement requirement of the
specific DSP implementation. 4 Other processors will/will not require
different input rearrangements which will affect results.

###  IV-B Transpose

Referring to the input DW vector as a virtual 2D matrix of
$\left[\big{\lceil}\frac{N}{4P}\big{\rceil}\right]\times[4P]$, the vector is
accessed with constant strides defined by the number of columns $[4P]$ as seen
in Fig. 4. This virtual 2D matrix is transposed such that indexes $j+4P$ are
the successors of indexes $j$. The CEVA-XC4500 DSP best transposes blocks of
$4\times 16[DWs]$, therefore a block of 4 rows, each of $16[DWs]$, is loaded
from the memory. Fig. 5 shows an example first block of the virtual 2D matrix
illustrated in Fig. 4. The $i^{th}$ block to be transposed is defined by four
rows of 16 consecutive DWs. The $k^{th}$ row of block $i$ starts at:

$16i+k\cdot 4P;k=0,1,2,3,1\leq i\leq\bigg{\lceil}\frac{4P}{16}\bigg{\rceil}$

and $\big{\lceil}\frac{4P}{16}\big{\rceil}$ is the number of blocks to be
transposed. The number of stuffed rows in the last block is
$16\cdot\big{\lceil}\frac{4P}{16}\big{\rceil}-4P$.

Figure 4: Example of a virtual 2D matrix for: $N=776;P=39$. Figure 5: The
first block of the virtual 2D matrix (highlighted in Fig. 4) before
transposition.

Fig. 6 shows the same block illustrated in Fig. 5. Each row is stored in a
continuous manner. The next row is stored at the end of the previous row. The
next block rows are stored at the end of the previous block.

Figure 6: Block in Fig. 5 after transposition.

Notes:

  * •

An important part of this phase is swapping the words of odd-indexed
($1,3,5,...,N-1$) DWs, meaning that in every odd-indexed DW the high-word and
the low-word are swapped. This swapping implements the permutation inside the
bit-couplet, as required. The CEVA-XC4500 DSP can swap word while loading any
DW, meaning the bit-level swapping is implemented without influencing
processing time or memory access.

  * •

The chosen vector size and parameter sets used in simulations require
transposing matrices of up-to eight rows. An expansion of this algorithm for
larger sized matrices might be possible but was not simulated.

###  IV-C Concatenation

This phase creates four vectors by concatenating rows of the transposed 2D
matrix of the previous phase. The four initial rows are chosen by four pre-
calculated offsets determined by the DVB-RCS2 standard parameters
$Q_{0},Q_{1},Q_{2},Q_{3}$ and $P$:

  * •

$\mathrm{offset}_{1}=3$

  * •

$\mathrm{offset}_{2}=(3+4Q_{1}+P)\ (\mathrm{mod}\ N)$

  * •

$\mathrm{offset}_{3}=(3+4(Q_{0}\cdot P+Q_{2})\ (\mathrm{mod}\ N)+2P)$

  * •

$\mathrm{offset}_{4}=(3+4(Q_{0}\cdot P+Q_{3})\ (\mathrm{mod}\ N)+3P)$

The strides taken from one row to the next are determined by
$G=mod_{N}(4P\cdot\big{\lceil}\frac{N}{4P}\big{\rceil})$. We continue and
construct each vector concatenating the row G from it. Knowing there are $4P$
rows in the transposed 2D matrix and 4 vectors, consequently there are $P-1$
concatenations to be done:

Figure 7: Concatenated $\mathrm{Vector}_{i};1\leq i\leq 4$.

Note: In Fig. 7 all row numbers must wrap around, keeping the index in range:
$0\leq row_{idx}\leq 4\cdot(P-1)$. Some parameter sets dictate a starting
point for concatenation that is not the first DW of the first chosen row. The
first DW for each $\mathrm{vector}_{i}$ is defined by:

$L_{k}=\bigg{\lceil}\frac{\mathrm{offset}_{k}}{4P}\bigg{\rceil};k=1,2,3,4$

Figure 8: Transposed 2D matrix (in Fig. 4) with offsets

Fig. 8 depicts the virtual 2D matrix illustrated in Fig. 4. In this example
$N=776;P=39;Q_{0}=7;Q_{1}=0;Q_{2}=0;Q_{3}=0$. Hence, using the expressions
above one may obtain: $\mathrm{offset}_{1}=3$; $\mathrm{offset}_{2}=42$;
$\mathrm{offset}_{3}=397$; $\mathrm{offset}_{4}=436$; $G=4$; $L_{1}=0$ (always
true since $\mathrm{offset}_{1}\equiv 3$); $L_{2}=0$, $L_{3}=2$; $L_{4}=2$.

Figure 9: Four concatenated vectors (generated from Fig. 8).

Fig. 9 shows an example of four result vectors in positions
$\mathrm{offset}_{1}-\mathrm{offset}_{4}$ which are the output of the
concatenation phase and serve as the input for the ordering phase.

###  IV-D Ordering

The ordering phase organises it’s input as one long vector. This is done by
loading the four vectors and transposing them using a similar transposition
process to the one described in the transpose phase. Here, every iteration
transposes blocks of $4\times 16$ [DWs] with the following differences:

  1. 1.

The rows transposed are now stored in an orderly fashion which match the
desired output vector (no need for reconstruction of the transposed matrix).

  2. 2.

The vectors are of same size $(\frac{N}{4})$, meaning that stuffed data will
exist only to fill the last block transposed.

Figure 10: Four concatenated vectors and matching output vector for
$N=776;P=39$.

Therefore, the number of blocks transposed is $\big{\lceil}\frac{N}{4\cdot
16}\big{\rceil}$, and the number of stuffed columns in the last block is
$16\cdot\bigg{\lceil}\frac{N}{4\cdot 16}\bigg{\rceil}-\frac{N}{4}$. Fig. 10
illustrates the $4\times\frac{N}{4}$ [DWs] matrix before reordering and the
output of this phase, $N$ [DWs] after permutation (example).

###  IV-E Word To Bit

The final operation of the algorithm is un-padding the output vector of the
ordering phase from the redundant zeros that were added to the original data
in the ”bit to word” phase.

This concludes the vector-wise permutation stage implementation. Given $N$ bit
couplets and matching parameters, following these phases, the result vector
will contain $N$ bit couplets after permutation. The described are compatible
with the DVB-RCS2 standard and are as generic as possible, having little data
adjustments to best-fit CEVA’s DSP can be reconstructed to fit any VP with
very few and minor changes to the overall algorithm.

The asymptotic tight bound computational complexity of the vector-wise
algorithm can be written as:

$O(N+\big{\lceil}\frac{N}{P}\big{\rceil}\cdot P)\simeq O(N)$

##  V Simulations And Results

Performance analysis was carried out using the CEVATOOLBOX profiler designed
for CEVA’s DSPs. The profiler creates a table detailing (1) the cycle count of
each function in the code and total cycle count results and (2) the number of
Read/Write operations performed in each function and total Read/Write results.
Unlike the approach of enhancing each function of a code/algorithm
individually, the proposed vector implementation revised the whole permutation
process. Hence, serial and vector algorithms are fundamentally incompatible,
therefore we consider only the overall results of the permutation stage
detailed in Tables II and sec.III.

The simulations were executed for 11 vector sizes and matching parameter sets
given in Table I. for each parameter set, we compared the serial and vector-
wise results and calculated the speedup using:

$Speedup=\frac{L_{old}-L_{new}}{L_{new}}\cdot 100\%.$

where $L_{old}$ and $L_{new}$ are the counts of Read, Write and Total Cycles
of the serial and vector-wise implementations, respectively.

TABLE II: Read/Write Profiler Results

Vector Size | Parallel | Speedup
---|---|---
N[bits] | Read | Write | Read | Write
56 | 761 | 666 | 11.7 | 109.6
152 | 621 | 570 | 82.1 | 158.6
236 | 739 | 802 | 200.0 | 143.8
384 | 822 | 719 | 284.7 | 230.4
432 | 988 | 1042 | 383.8 | 199.0
492 | 1279 | 997 | 623.4 | 409.5
520 | 998 | 1042 | 438.5 | 222.2
776 | 1217 | 928 | 589.9 | 406.4
1056 | 1487 | 1109 | 307.1 | 229.8
1192 | 1487 | 1109 | 734.5 | 484.4
2396 | 1612 | 1189 | 295.2 | 219.6

TABLE III: Total Cycles Profiler Results

Vector Size | Parallel | Speedup
---|---|---
N[bits] | Total Cycles | Total Cycles
56 | 172239 | 6.000
152 | 132264 | 63.70
236 | 175387 | 87.70
384 | 178026 | 141.7
432 | 238599 | 153.1
492 | 270537 | 296.2
520 | 238599 | 177.8
776 | 255186 | 284.1
1056 | 313108 | 134.9
1192 | 313108 | 339.9
2396 | 339651 | 124.7

The results of Tables II and III are also presented in Fig. 11.

Figure 11: Speedup results of Read, Write and Total Cycles by Vector Size

We notice a positive speed-up factor for all vector sizes. The reason for the
non-monotonic plots is the different execution of the vector algorithm caused
by the given parameters of the each vector size and not only by the vector
size itself. Re-analyzing the vector algorithm phases by the parameters we see
that for example the number of blocks transposed is a factor of P and can
shorten/lengthen the runtime of the application by reducing/adding blocks to
transpose. See Sec. VI for more.

##  VI Conclusion and Future Research

###  VI-A Conclusion

The vector-wise algorithm implementation out performs the serial
implementation for all tested cases. Speedups were achieved even though the
DSP hardware required padding and unpadding of the original data and final
result respectively. A hardware better fit for the parallel algorithm,
enabling the Transpose, Concatenation and Ordering phases to be executed on
bits rather than words will improve the superiority of the parallel algorithm
over the serial implementation.

###  VI-B Future Research

Research possibilities include different paths:

  1. 1.

Additional permutation algorithms can and be tested on the DSP.

  2. 2.

The parallel algorithm can be implemented or simulated on other processors
(e.g. GPUs, CPUs, etc.).

  3. 3.

The encoding performed by this application is targeted at achieving maximum
channel capacity with minimum bit error rate (BER). The parameters used in
this research (shown in Table I) were given as constants. A wider approach
would be to optimize vector length with matching parameters, while considering
the algorithm presented and it’s vector properties, all this while striving to
achieve better BER and throughput.

## References

  * [1] C. E. Shannon, “A mathematical theory of communication,” _The Bell System Technical Journal_ , vol. 27, no. 3, pp. 379–423, July 1948.
  * [2] C. Berrou, A. Glavieux, and P. Thitimajshima, “Near shannon limit error-correcting coding and decoding: Turbo-codes. 1,” in _Communications, 1993. ICC ’93 Geneva. Technical Program, Conference Record, IEEE International Conference on_ , vol. 2, May 1993, pp. 1064–1070 vol.2.
  * [3] Aliazam-Abbasfar, _Turbo-like Codes, Design for High Speed Decoding_ , 1st ed., ser. 1. Springer, 1 2007, vol. 1.
  * [4] C. Douillard, M. Jezequel, C. Berrou, J. Tousch, N. Pham, and N. Brengarth, “The turbo code standard for dvb-rcs,” in _Proc. of the International Symposium on Turbo Codes and Related Topics (ISTC)_ , 2000, pp. 535–538.
  * [5] G. Oletu and P. Rapajic, “The performance of turbo codes for wireless communication systems,” in _Computer Research and Development (ICCRD), 2011 3rd International Conference on_ , vol. 4, March 2011, pp. 346–349.
  * [6] T. Ilnseher, F. Kienle, C. Weis, and N. Wehn, “A 2.15gbit/s turbo code decoder for lte advanced base station applications,” in _Turbo Codes and Iterative Information Processing (ISTC), 2012 7th International Symposium on_ , Aug 2012, pp. 21–25.
  * [7] P. Bahirgonde and S. Dixit, “Ber analysis of turbo code interleaver,” in _International Journal of Computer Applications Volume 126 - No.14_ , Sep 2015\.
  * [8] J. P. Woodard and L. Hanzo, “Comparative study of turbo decoding techniques: an overview,” _IEEE Transactions on Vehicular Technology_ , vol. 49, no. 6, pp. 2208–2233, Nov 2000.
  * [9] O. Nordmark, “Turbo code performance analysis using hardware acceleration,” 2016\.
  * [10] A. Kudriavtsev and P. Kogge, “Generation of permutations for simd processors,” in _ACM SIGPLAN Notices_ , vol. 40, no. 7. ACM, 2005, pp. 147–156.
  * [11] R. Barik, J. Zhao, and V. Sarkar, “Efficient selection of vector instructions using dynamic programming,” in _Microarchitecture (MICRO), 2010 43rd Annual IEEE/ACM International Symposium on_. IEEE, 2010, pp. 201–212.


# Approximation and bounds for the Wallis ratio

Xu You

###### Abstract

In this paper, we present an improved continued fraction approximation of the
Wallis ratio. This approximation is fast in comparison with the recently
discovered asymptotic series. We also establish the double-side inequality
related to this approximation. Finally, some numerical computations are
provided for demonstrating the superiority of our approximation.

0002010 Mathematics Subject Classification: 33B15, 26A48, 26D07000Key words
and phrases: Wallis ratio, Gamma function, Inequalities, Multiple-correction
method

##  1 Introduction

The Wallis ratio is defined as

$W(n)=\frac{(2n-1)!!}{(2n)!!}=\frac{1}{\sqrt{\pi}}\frac{\Gamma(n+\frac{1}{2})}{\Gamma(n+1)},$

where $\Gamma$ is the classical Euler gamma function which may be defined by

$\Gamma(x)=\int_{0}^{\infty}{t^{x-1}e^{-t}}dt,\textmd{Re}(x)>0.$

The study and applications of $W(n)$ have a long history, a large amount of
literature, and a lot of new results. For detailed information, please refer
these papers [1, 2, 3, 4] and references cited therein.

Chen and Qi [1] presented the following inequalities for the Wallis ratio for
every natural number $n$:

(1.1)
$\displaystyle\frac{1}{\sqrt{\pi(n+\frac{4}{\pi}-1)}}\leq\frac{(2n-1)!!}{(2n)!!}\leq\frac{1}{\sqrt{\pi(n+\frac{1}{4})}},$

where the constants $\frac{4}{\pi}-1$ and $\frac{1}{4}$ are the best possible.

Guo, Xu and Qi proved in [6] that the double inequality

(1.2)
$\displaystyle\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2n}\right)^{n}\frac{\sqrt{n-1}}{n}<W(n)\leq\frac{4}{3}\left(1-\frac{1}{2n}\right)^{n}\frac{\sqrt{n-1}}{n},$

for $n\geq 2$ is valid and sharp in the sense that the constants
$\sqrt{\frac{e}{\pi}}$ and $\frac{4}{3}$ are best possible. They also proposed
the approximation formula

(1.3) $\displaystyle
W(n)\sim\chi(n):=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2n}\right)^{n}\frac{\sqrt{n-1}}{n},n\rightarrow\infty.$

Recently, Qi and Mortici [5] improved the approximation formula (1.3) as
following,

(1.4) $\displaystyle
W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}.$

Motivated by these works, in this paper we will apply the _multiple-correction
method_ [7, 8, 9] to construct an improved continued fraction asymptotic
expansion for the Wallis ratio as follows:

######  Theorem 1.

For the Wallis ratio $W(n)=\frac{(2n-1)!!}{(2n)!!}$, we have

(1.5) $\displaystyle
W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}+\frac{a_{2}}{n+b_{2}+\frac{a_{3}}{n+b_{3}}}}\right),$

where
$a_{1}=\frac{1}{144},b_{1}=\frac{1}{60};a_{2}=\frac{781}{3600},b_{2}=-\frac{4309}{109340};a_{3}=\frac{51396085}{89664267},b_{3}=\frac{25682346121}{449571834712}.$

Using Theorem 1, we provide some inequalities for the Wallis ratio.

######  Theorem 2.

For every integer $n>1$, it holds:

(1.6)
$\displaystyle\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}+\frac{a_{2}}{n+b_{2}+\frac{a_{3}}{n+b_{3}}}}\right)$
$\displaystyle>W(n)>\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}+\frac{a_{2}}{n+b_{2}}}\right),$

where
$a_{1}=\frac{1}{144},b_{1}=\frac{1}{60};a_{2}=\frac{781}{3600},b_{2}=-\frac{4309}{109340};a_{3}=\frac{51396085}{89664267},b_{3}=\frac{25682346121}{449571834712}.$

To obtain Theorem 1, we need the following lemma which was used in [10, 11,
12] and is very useful for constructing asymptotic expansions.

######  Lemma 1.

If the sequence $(x_{n})_{n\in\mathbb{N}}$ is convergent to zero and there
exists the limit

(1.7)
$\displaystyle\lim_{n\rightarrow+\infty}n^{s}(x_{n}-x_{n+1})=l\in[-\infty,+\infty]$

with $s>1$, then

(1.8) $\displaystyle\lim_{n\rightarrow+\infty}n^{s-1}x_{n}=\frac{l}{s-1}.$

Lemma 1 was proved by Mortici in [10]. From Lemma 1, we can see that the speed
of convergence of the sequences $(x_{n})_{n\in\mathbb{N}}$ increases together
with the values $s$ satisfying (1.8).

The rest of this paper is arranged as follows. In section 2, we will apply the
_multiple-correction method_ to construct a new asymptotic expansion for the
Wallis ratio and prove Theorem 1 by the _multiple-correction method_. In
section 3, we established the double-side inequality for the Wallis ratio. In
section 4, we give some numerical computations which demonstrate the
superiority of our new series over some formulas found recently.

##  2 Proof of Theorem 1

According to the argument of Theorem 5.1 in [5] , we can introduce a sequence
$(u(n))_{n\geq 1}$ by the relation

(2.1) $\displaystyle
W(n)=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp
u(n),$

and to say that an approximation
$W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}$
is better if the speed of convergence of $u(n)$ is higher.

(Step 1) The initial-correction. When $n\rightarrow\infty$, we define a
sequence $(u_{0}(n))_{n\geq 1}$

(2.2) $\displaystyle
W(n)=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp
u_{0}(n).$

From (2.2), we have

(2.3) $\displaystyle u_{0}(n)=\ln
W(n)-\ln\sqrt{\frac{e}{\pi}}-\left(n+\frac{1}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)-\frac{1}{2}\ln\frac{1}{n}.$

Thus,

(2.4) $\displaystyle u_{0}(n)-u_{0}(n+1)=$
$\displaystyle\ln\frac{2n+2}{2n+1}-\left(n+\frac{1}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)-\frac{1}{2}\ln\frac{1}{n}$
$\displaystyle+\left(n+\frac{4}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{4}{3}\right)}\right)+\frac{1}{2}\ln\frac{1}{n+1}.$

Developing (2.4) into power series expansion in $1/n$, we have

(2.5) $\displaystyle
u_{0}(n)-u_{0}(n+1)=\frac{1}{48}\frac{1}{n^{4}}+O(\frac{1}{n^{5}}),$

By Lemma 1, we know that the rate of convergence of the sequence
$(u_{0}(n))_{n\geq 1}$ is $n^{-3}$.

(Step 2) The first-correction. We define the sequence $(u_{1}(n))_{n\geq 1}$
by the relation

(2.6) $\displaystyle
W(n)=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}}\right)\exp
u_{1}(n).$

From (2.6), we have

(2.7) $\displaystyle u_{1}(n)-u_{1}(n+1)=$
$\displaystyle\ln\frac{2n+2}{2n+1}-\left(n+\frac{1}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)-\frac{1}{2}\ln\frac{1}{n}+\frac{1}{2}\ln\frac{1}{n+1}$
$\displaystyle+\left(n+\frac{4}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{4}{3}\right)}\right)-\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}}+\frac{1}{(n+1)^{2}}\frac{a_{1}}{n+1+b_{1}}.$

Developing (2.7) into power series expansion in $1/n$, we have

(2.8) $\displaystyle
u_{1}(n)-u_{1}(n+1)=\left(\frac{1}{48}-3a_{1}\right)\frac{1}{n^{4}}+\left(-\frac{91}{2160}+a_{1}(6+4b_{1})\right)\frac{1}{n^{5}}+O(\frac{1}{n^{6}}),$

By Lemma 1, we know that the fastest possible sequence $(u_{1}(n))_{n\geq 1}$
is obtained as the first item on the right of (2.8) vanishes. So taking
$a_{1}=\frac{1}{144},b_{1}=\frac{1}{60}$, we can get the rate of convergence
of the sequence $(u_{1}(n))_{n\geq 1}$ is at least $n^{-5}$.

(Step 3) The second-correction. We define the sequence $(u_{2}(n))_{n\geq 1}$
by the relation

(2.9) $\displaystyle
W(n)=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}+\frac{a_{2}}{n+b_{2}}}\right)\exp
u_{2}(n).$

Using the same method as above, we obtain that the sequence $(u_{2}(n))_{n\geq
1}$ converges fasted only if
$a_{2}=\frac{781}{3600},b_{2}=-\frac{4309}{109340}$.

(Step 4) The third-correction. Similarly, define the sequence
$(u_{3}(n))_{n\geq 1}$ by the relation

(2.10) $\displaystyle
W(n)=\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{a_{1}}{n+b_{1}+\frac{a_{2}}{n+b_{2}+\frac{a_{3}}{n+b_{3}}}}\right)\exp
u_{3}(n)$

Using the same method as above, we obtain that the sequence $(u_{3}(n))_{n\geq
1}$ converges fasted only if
$a_{3}=\frac{51396085}{89664267},b3=\frac{25682346121}{449571834712}$.

The new asymptotic (1.5) is obtained.

##  3 Proof of Theorem 2

The double-side inequality (1.6) may be written as follows:

$f(n)=\ln
W(n)-\frac{1}{2}+\frac{1}{2}\ln\pi-\left(n+\frac{1}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)-\frac{1}{2}\ln\frac{1}{n}-\frac{1}{n^{2}}\frac{\frac{1}{144}}{n+\frac{1}{60}+\frac{\frac{781}{3600}}{n-\frac{4309}{109340}}}$

and

$\displaystyle g(n)=$ $\displaystyle\ln
W(n)-\frac{1}{2}+\frac{1}{2}\ln\pi-\left(n+\frac{1}{3}\right)\ln\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)$
$\displaystyle-\frac{1}{2}\ln\frac{1}{n}-\frac{1}{n^{2}}\frac{\frac{1}{144}}{n+\frac{1}{60}+\frac{\frac{781}{3600}}{n-\frac{4309}{109340}+\frac{\frac{51396085}{89664267}}{n+\frac{25682346121}{449571834712}}}}.$

Suppose $F(n)=f(n+1)-f(n)$ and $G(n)=g(n+1)-g(n)$. For every $x>1$, we can get

(3.1) $\displaystyle
F^{\prime\prime}(x)=\frac{A(x-1)}{160x^{4}(1+x)^{4}(1+2x)^{2}(1+3x)(4+3x)(6x-1)^{2}(5+6x)^{2}\Psi_{1}^{3}(x;2)\Psi_{2}^{3}(x;2)}<0$

and

(3.2) $\displaystyle
G^{\prime\prime}(x)=\frac{B(x)}{32x^{4}(1+x)^{4}(1+2x)^{2}(1+3x)(4+3x)(-1+6x)^{2}(5+6x)^{2}\Psi_{1}^{3}(x;3)\Psi_{2}^{3}(x;3)}>0,$

where

$\displaystyle\Psi_{1}(x;2)=$ $\displaystyle 10642-1119x+49203x^{2},$
$\displaystyle\Psi_{2}(x;2)=$ $\displaystyle 58726+97287x+49203x^{2},$
$\displaystyle\Psi_{1}(x;3)=$ $\displaystyle
113505180+4083418255x+178132605x^{2}+5180725368x^{3}$
$\displaystyle\Psi_{2}(x;3)=$ $\displaystyle
9555781408+19981859569x+15720308709x^{2}+5180725368x^{3},$

$A(x)=-461003,...,512n^{18}-...$ is a polynomial of 18th degree with all
negative coefficients and $B(x)=212876,...,696x^{22}+...$ is a polynomial of
22th degree with all positive coefficients.

It show that $F(x)$ is strictly concave and $G(x)$ is strictly convex on
$(0,\infty)$. According to Theorem 1, when $n\rightarrow\infty$, it holds that
$\lim_{n\rightarrow\infty}f(n)=\lim_{n\rightarrow\infty}g(n)=0$; thus
$\lim_{n\rightarrow\infty}F(n)=\lim_{n\rightarrow\infty}G(n)=0$. As a result,
we can make sure that $F(x)<0$ and $G(x)>0$ on $(0,\infty)$. Consequently, the
sequence $f(n)$ is strictly increasing and $g(n)$ is strictly decreasing while
they both converge to $0$. As a result, we conclude that $f(n)>0$, and
$g(n)<0$ for every integer $n>1$.

The proof of Theorem 2 is complete.

##  4 Numerical computations

In this section, we give Table 1 to demonstrate the superiority of our new
series respectively. From what has been discussed above, we found out the new
asymptotic function as follows:

(4.1) $\displaystyle
W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\frac{1}{\sqrt{n}}\exp\left(\frac{1}{n^{2}}\frac{\frac{1}{144}}{n+\frac{1}{60}+\frac{\frac{781}{3600}}{n-\frac{4309}{109340}+\frac{\frac{51396085}{89664267}}{n+\frac{25682346121}{449571834712}}}}\right)=\sigma(n).$

Chen and Qi [2] gave:

(4.2) $\displaystyle
W(n)\sim\frac{1}{\sqrt{\pi\left(n+\frac{1}{4}\right)}}=\alpha(n).$

Qi and Mortici [5] gave the improved formula:

(4.3) $\displaystyle
W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2n}\right)^{n}\sqrt{\frac{1}{n}}\exp\left(\frac{1}{24n^{2}}+\frac{1}{48n^{3}}+\frac{1}{160n^{4}}+\frac{1}{960n^{5}}\right)=\beta(n)$

and

(4.4) $\displaystyle
W(n)\sim\sqrt{\frac{e}{\pi}}\left(1-\frac{1}{2\left(n+\frac{1}{3}\right)}\right)^{n+\frac{1}{3}}\sqrt{\frac{1}{n}}=\gamma(n).$

We can easily observe that the new formula converges fastest of the other
three formulas.

Table 1: Simulations for $\alpha(n),\beta(n),\gamma(n)$ and $\sigma(n)$ n | $\frac{W(n)-\alpha(n)}{W(n)}$ | $\frac{W(n)-\beta(n)}{W(n)}$ | $\frac{W(n)-\gamma(n)}{W(n)}$ | $\frac{W(n)-\sigma(n)}{W(n)}$
---|---|---|---|---
50 | $-6.1876\times 10^{-6}$ | $7.3576\times 10^{-14}$ | $5.5532\times 10^{-8}$ | $-3.8082\times 10^{-19}$
500 | $-6.2438\times 10^{-8}$ | $7.1643\times 10^{-20}$ | $5.5554\times 10^{-11}$ | $-3.8138\times 10^{-28}$
1000 | $-1.5617\times 10^{-8}$ | $1.1177\times 10^{-21}$ | $6.9443\times 10^{-12}$ | $-7.4489\times 10^{-31}$
2000 | $-3.9053\times 10^{-9}$ | $1.7452\times 10^{-23}$ | $8.6805\times 10^{-13}$ | $-1.4549\times 10^{-33}$

##  5 Acknowledgements

This work was supported by the National Natural Science Foundation of China
(Grant No. 61403034), Beijing Municipal Commission of Education Science and
Technology Program (KM201510017002).

## References

  * [1] C. P. Chen, F. Qi, Completely monotonic function associated with the gamma functions and proof of Wallis inequality, Tamkang J. Math. 36(4) (2005) 303-307.
  * [2] C. P. Chen, F. Qi, The best bounds in Wallis inequality, Proc. Amer. Math. Soc. 133 (2) (2005) 397 C401.
  * [3] F. Qi, Integral representations and complete monotonicity related to the remainder of Burnside s formula for the gamma function, J. Comput. Appl. Math. 268 (2014) 155 C167.
  * [4] F. Qi, L. H. Cui, S.L. Xu, Some inequalities constructed by Tchebysheff s integral inequality, Math. Inequal. Appl. 2 (4) (1999) 517 C528.
  * [5] F. Qi, C. Mortici, Some best approximation formulas and inequalities for the Wallis ratio, Appl. Math. Comput. 253 (2015) 363-368.
  * [6] S. Guo, J. G. Xu, F. Qi, Some exact constants for the approximation of the quantity in the Wallis formula, J. Inequal. Appl. 2013 (67) (2013). 7 pages.
  * [7] X. D. Cao, H. M. Xu, X. You, Multiple-correction and faster approximation. J. Number Theory. 149 (2015), 327-350.
  * [8] X. D. Cao, Multiple-correction and continued fraction approximation J. Math. Anal. Appl. 424 (2015), 1425-1446.
  * [9] X. D. Cao, X. You, Multiple-correction and continued fraction approximation(II). Appl. Math. Comput. 261 (2015) 192-205.
  * [10] C. Mortici, Product approximations via asymptotic integration, Amer. Math. Monthly 117(5) (2010) 434-441.
  * [11] C. Mortici, New improvements of the Stirling formula, Appl. Math. Comput. 217(2) (2010) 699-704.
  * [12] C. Mortici, A continued fraction approximation of the gamma function, J. Math. Anal. Appl. 402 (2013) 405-410.

Xu You

1\. Department of Mathematics and Physics, Beijing Institute of Petrochemical
Technology,

Beijing 102617, P. R. China

e-mail<EMAIL_ADDRESS>