text:
# A Spectrophotometric Method to Determine the Inclination of Class I Objects

## 1 Introduction

Low mass pre-main-sequence stars are classified into four categories, Class 0, I, II, and III, mainly according to the shape of their spectral energy distributions (SEDs) (Lada & Wilking 1984; Adams, Lada, & Shu 1987, André, Ward-Thompson, & Barsony 1993). The shape of SED is strongly dependent on the structure surrounding a central star, such as an envelope, an outflow, a disk, and so on. For example, since Class 0 and Class I objects show a peak in their SEDs in a sub-millimeter range, they are thought to be surrounded by a cool massive envelope, which heavily obscures a hot central part of the object. It is considered generally that the order of the classification, 0, I, II, and III corresponds to an evolutional sequence of a low mass pre-main-sequence star. Class 0 and I objects are thought to be protostars, which are enshrouded by a thick envelope, and to be younger than T Tauri stars, which are classified as Class II and III objects.

A protostar is not a spherically symmetric system as implied by the presence of a bipolar outflow. This asymmetry could come from a structure such as a central star with a disk system, which seems to be the progenitor of a planetary system. Thus, a careful scrutiny of the structure of a protostar system is important to understand the formation process of a star and/or planetary system. However, it is not easy to reveal the structure directly by imaging observations, because optical imaging is significantly subject to the optically thick envelope and also radio imaging is strongly limited by the resolution of the technique. In contrast, a spectrophotometric method using SEDs seems to be effective to investigate a large number of protostar systems (Kenyon, Calvet, & Hartmann 1993; Men’shchikov, Henning, & Fischer 1999; Nakamoto & Kikuchi 1999). For instance, Kenyon et al. (1993) modeled a protostar with incorporating spherically symmetric radiative transfer and simulated its SED. They fitted simulated SEDs with observed SEDs and inferred the structure of Class I objects. The accuracy of the spectrophotometric method depends on the accuracy of the model of protostar, which is characterized by the density and temperature distributions in the envelope, the inclination of the system, and so forth. In particular, the inclination of the system is an important piece of information to probe the inner structure of the protostar. It can affect the SED significantly through 2D behavior of radiative transfer. Hence, the inclination is of essential importance in a realistic model of protostars. In other words, there may be a possibility to estimate the inclination using the SED. So far, the inclination of a protostar has been estimated by the geometrical feature of molecular outflow lobes. However, the geometrical method requires high spatial resolution of observations and long integration time. Thus, it is not effective for a statistical study. On the other hand, if we can estimate the inclination of the protostar solely by the SED, it would be quite effective to investigate the protostar structure.

In this paper, we propose a new method to determine the inclination of Class I objects using their SEDs. We solve the radiative equilibrium for the 2D axisymmetric density distribution with incorporating 2D radiative transfer and simulate SEDs for the 2D protostar model. Then, it is shown that the ratio of the quantity $`L_{\mathrm{SED}}`$, which is the integration of the SED over frequencies, to the peak flux of the SED, $`(\nu L_\nu )_{\mathrm{max}}`$, can be a good measure of the inclination. The accuracy and the applicability of the present method is discussed as well.

This paper is organized as follows. In section 2 and 3, our model of a protostar and calculation procedure of the new method will be described. The new method will be introduced in section 4. The validity, reliability, and accuracy of the new method is also discussed in this section. The conclusions of this paper will be given in section 5.

## 2 Model

The present two-dimensional axisymmetric model for a Class I object consists of three components; a central star, a circumstellar disk, and an envelope surrounding the former two components. The main opacity source is dust grains. It is assumed that the dust grains are well mixed with hydrogen molecule gas and the mass ratio to the hydrogen gas is $`0.01`$. The opacity of dust grains as a function of the frequency is taken from Miyake & Nakagawa (1993). Physical parameters that characterize our model are summarized in Table 1.

The central star is only the energy source and radiative heating is considered. Compressional heating, viscous heating, chemical heating, and any kind of heating sources other than the radiation from the central star are ignored, since they are negligibly small in almost all the regions of interest. The central star is assumed to have the surface temperature $`T_{}`$ and emit the blackbody radiation with that temperature. The luminosity of the central star is given by $`L_{}`$. The mass of the central star is denoted by $`M_{}`$.

Our model contains a circumstellar disk, although the presence of the disk has not been confirmed observationally in Class I objects. Its inner and outer radii are assumed to be $`0.1\mathrm{AU}`$ and $`100\mathrm{A}\mathrm{U}`$, respectively. The surface density distribution of the disk is given by $`\mathrm{\Sigma }(r)=\mathrm{\Sigma }_1(r/1\mathrm{A}\mathrm{U})^p`$, where $`r`$ is the radius in cylindrical coordinates from the central star, $`\mathrm{\Sigma }_1`$ is surface density at $`1\mathrm{A}\mathrm{U}`$, and $`p`$ is the index of the power, respectively. We adopt $`p=1.5`$ which is a standard value for the minimum mass solar nebula model (Hayashi, Nakazawa, & Nakagawa 1985). The value of $`p`$ does not affect the SEDs significantly. Vertical structure in the disk is determined by the hydrostatic equilibrium between gravity by the central star and thermal pressure.

We assume the outer radius of the envelope to be $`1,000\mathrm{A}\mathrm{U}`$. Contributions to SED from the outside of that region are neglected. This assumption seems appropriate, since the observed region of a Class I object is as small as a few thousand AU. In the envelope, the density distribution is given by

$$\rho (R)=\rho _1\left(\frac{R}{1\mathrm{A}\mathrm{U}}\right)^q,\text{ for }1\mathrm{A}\mathrm{U}<R<1,000\mathrm{A}\mathrm{U}\text{}$$
(1)
where $`R`$ is distance from the central star in spherical coordinates, $`\rho _1`$ is mass density at $`1\mathrm{A}\mathrm{U}`$, and $`q`$ is the index of the power, respectively. Here, we adopt $`q=1.5`$ because it is the value expected for the free-falling envelope (Shu 1977). Observationally, it is known that Class I objects usually have bipolar outflows, and it is suggested that the outflow excavates the envelope and affects the SED considerably. In order to take into account the effects of the bipolar outflow, we situate a cavity along the symmetric axis. The cavity is assumed to be a conical shape with a half opening angle $`\theta _{\mathrm{bp}}`$, which is the angle between the symmetric axis and the boundary of the cavity. The density in the cavity is assumed to be 0.01 times smaller than the value which is given by eq.(1). Thus, the density in the envelope is given by

$$\rho (R,\theta )=\{\begin{array}{cc}\rho (R)\hfill & \text{ for }\theta >\theta _{\mathrm{bp}}\text{}\hfill \\ 0.01\times \rho (R)\hfill & \text{ for }\theta <\theta _{\mathrm{bp}}\text{}\hfill \end{array}$$
(2)
## 3 2D Radiative Equilibrium Calculations

### 3.1 Numerical Method

In the first step, we obtain the temperature distribution in the model under the condition of the radiative equilibrium; $`_0^{\mathrm{}}\chi _\nu ^{\mathrm{abs}}J_\nu 𝑑\nu =_0^{\mathrm{}}\chi _\nu ^{\mathrm{abs}}B_\nu 𝑑\nu `$, where $`\chi _\nu ^{\mathrm{abs}}`$ is the absorption coefficient, $`J_\nu `$ is the mean intensity, and $`B_\nu `$ is the Planck function, respectively. To calculate the radiative equilibrium, we solve 2D radiative transfer by a numerical scheme based on the so-called Variable Eddington Factor method (Stone, Mihalas, & Norman 1992; Kikuchi, Nakamoto, & Ogochi 2002). The primary difference of our algorithm from that by Stone et al. is that we treat a frequency dependence of radiative transfer accurately. In our method, the radiation field is separated into two parts; one is the direct component from the central star and the other is the diffuse component. The direct component is calculated straightforwardly. The diffuse component is obtained by solving moment equations (Stone et al. 1992), which are extended to a non-gray version in our method, coupled with the energy equation. To close the moment equations, the variable Eddington factor is introduced, which is determined by solving the radiative transfer equation taking the scattering effect into account. We expand the scattering phase function by the Legendre function to the 2nd order. The set of equations are integrated with time. When the temporal change of the temperature distribution in the calculation region becomes small enough, we regard it as an equilibrium state. In the second step, the SED of the object is simulated by the direct ray-tracing with the radiative equilibrium temperature distribution.

Our numerical code is based on cylindrical coordinates. We use $`50\times 50`$ grids for $`r`$-$`z`$ coordinates in space, and $`50\times 100`$ grids for azimuthal and zenith angles to express directions in which the intensity propagates. Calculations with double number of grids in each dimension (16 times larger calculations) showed that the difference of results was within about 10% in $`f_\mathrm{L}`$, which will be defined in eq.(4). The frequency space is divided into 51 grids.

### 3.2 SED calculation

Simulated SEDs are illustrated in Fig. 1. It is seen that SEDs are sensitive to $`i`$, particularly in the frequency range higher than the peak frequency $`\nu _{\mathrm{peak}}`$. This change reflects the asymmetry of the density and temperature distribution. The optical depth along the line-of-sight from an observer to the central star, $`\tau _\mathrm{c}`$, is $`1`$ at $`\nu _{\mathrm{peak}}`$, while $`\tau _\mathrm{c}>1`$ for $`\nu >\nu _{\mathrm{peak}}`$ and $`\tau _\mathrm{c}<1`$ for $`\nu <\nu _{\mathrm{peak}}`$. Thus, for $`\nu <\nu _{\mathrm{peak}}`$, the dependence on the inclination is within a factor of $`23`$ and the emergent flux is the superposition of the thermal radiation from dust grains, which is determined by the temperature and the total amount of dust, regardless of density distributions. In contrast, for $`\nu >\nu _{\mathrm{peak}}`$, the anisotropy of flux is very large (2-4 orders of magnitude), and the emergent flux is dominated by the scattered radiation originating from the central star and inner disk, which is strongly affected by the density distribution of the envelope. Therefore, the flux in the pole-on ($`i=0^{}`$) direction is large, while the flux in the edge-on ($`i=90^{}`$) direction is small.

### 3.3 Comparison between Semi and Full 2D Calculations

To see the importance of full 2D radiative transfer calculations, we performed approximate semi 2D calculations, which were carried out by Kenyon et al. (1993). They derived the temperature distribution in the envelope from radiative equilibrium, but they assumed the spherically symmetric density distribution when they calculated the radiative equilibrium, so they obtained spherically symmetric temperature distribution. To do that, first the spherically averaged density profile from the 2D axially symmetric density distribution is obtained, and then the 1D spherically symmetric temperature profile is calculated. Using the spherical temperature distribution and the non-spherical 2D density distribution, the emergent SEDs from the object are obtained by ray-tracing . Following their procedure, we also calculated the SED from the semi 2D calculations and compared the results with our full 2D calculations.

In Fig. 2, we compare two temperature distributions derived from two procedures. The temperature distribution obtained by the full 2D calculations is not spherically symmetric at all. The most noticeable difference between two temperatures is found in the disk and the outflow regions. The temperature in the disk in the full 2D calculation is much lower than that of the spherical temperature, because the disk has higher density than the spherically averaged density. We also find that the temperature in a shaded region, which is behind the disk with respect to the central star, is lower than the spherical temperature, though the density in that region is not so different from spherically averaged one. This is because that the total amount of radiation from the central star to the region is reduced by the efficient disk absorption. In contrast, the temperature in the outflow region in the full 2D calculation is higher than spherical one, because the density in the region is lower than the spherically averaged one. Consequently, it is obvious that the temperature distribution obtained by the full 2D radiative equilibrium calculation is quite different from spherically symmetric distribution.

Emergent SEDs based on the temperature distributions displayed in Fig. 2 are shown in Fig. 3. For each inclination, it is seen that (1) the flux derived from the full 2D calculations in a frequency range from $`10^{12}\text{ Hz}`$ to $`10^{13}\text{ Hz}`$ is lower than that of the semi 2D calculations, and (2) in a range from $`10^{13}\text{ Hz}`$ to $`10^{15}\text{ Hz}`$, the full 2D calculation flux is higher than the semi 2D calculation one. The first difference is due to the decrease of the envelope temperature compared to the spherical temperature. Since the optical depth for the radiation at this frequency from the central star to the observer is around unity, the flux in this range is expected to be in proportion to the dust temperature in the envelope. On the other hand, the second difference is caused by the different mean intensity $`J_\nu `$ especially in the outflow region, because flux in this frequency range is mainly determined by the scattering, which is evaluated by $`\chi _\nu ^{\mathrm{sca}}J_\nu `$, where $`\chi _\nu ^{\mathrm{sca}}`$ is scattering coefficient, and is very anisotropic in the non-spherical density distribution cases. We can see from Figs. 2 and 3 that the full 2D radiative equilibrium calculations are indispensable to study the detailed structure of Class 0/I objects using their SEDs and/or images.

## 4 Estimation of Inclination Angle

### 4.1 A New Inclination Indicator: $`f_\mathrm{L}`$

The inclination dependence of SED provides a tool to infer the inclination itself. We introduce the emergent luminosity $`L_{\mathrm{SED}}`$, which is defined by

$$L_{\mathrm{SED}}=4\pi D^2_0^{\mathrm{}}(\nu F_\nu )d\mathrm{log}\nu ,$$
(3)
as a tracer for the change of SED, where $`D`$ is the distance to the target object and $`F_\nu `$ is the observed flux at frequency $`\nu `$. Also, using the peak flux $`(\nu L_\nu )_{\mathrm{max}}=(\nu 4\pi D^2F_\nu )_{\mathrm{max}}`$ (flux at $`\nu _{\mathrm{peak}}`$), we define a ratio $`f_\mathrm{L}`$ by

$$f_\mathrm{L}=\frac{L_{\mathrm{SED}}}{(\nu L_\nu )_{\mathrm{max}}}.$$
(4)
(Note that $`L_{\mathrm{SED}}`$, $`(\nu L_\nu )_{\mathrm{max}}`$, and $`f_\mathrm{L}`$ are all observable quantities. We can obtain them from observed SEDs.) We show the change of $`f_\mathrm{L}`$ against $`i`$ in Fig. 4. Fig. 4(a) indicates that $`i`$ can be determined if $`f_\mathrm{L}`$ is evaluated from observed SED for the angle of $`i\theta _{\mathrm{bp}}`$. But, when $`i`$ is in a range $`i<\theta _{\mathrm{bp}}`$, $`f_\mathrm{L}`$ is almost independent of $`i`$. This means that SED does not change if we observe a protostar through the bipolar cavity. For $`i90^{}`$, $`f_\mathrm{L}`$ is also independent of $`i`$. It corresponds to the case that we observe a protostar through the disk.

The level of $`f_\mathrm{L}`$ depends on $`\theta _{\mathrm{bp}}`$. It implies that $`\theta _{\mathrm{bp}}`$ is one of the important parameters to determine the anisotropy of a protostar system. The ratio $`f_\mathrm{L}`$ increases for all $`i`$ as $`\theta _{\mathrm{bp}}`$ increases, and it is worth noting that there exist maximum and minimum of $`f_\mathrm{L}`$ for each $`\theta _{\mathrm{bp}}`$.

We calculated $`f_\mathrm{L}`$ for full 2D and semi 2D results shown in Figs. 2 and 3 to clarify the importance of the full 2D calculation. These results are shown in Fig. 5. In the semi 2D calculation, the value of $`f_\mathrm{L}`$ is lower than that of the full 2D calculation, especially for the small inclination. It is because that the peak fluxes $`(\nu L_\nu )_{\mathrm{max}}`$ are different by a factor of $`2`$ while the difference of observed luminosities $`L_{\mathrm{SED}}`$ is so small. For the large inclination, the differences of $`(\nu L_\nu )_{\mathrm{max}}`$ and $`L_{\mathrm{SED}}`$ are canceled each other so that the values of $`f_\mathrm{L}`$ are almost the same.

### 4.2 Test of the Method

We test our new method with a real object TMC1A (IRAS 04365+2535). TMC1A is a Class I object (e.g. Kenyon et al. 1993; Chandler et al. 1998) and has an outflow with the half opening angle $`\theta _{\mathrm{bp}}20^{}`$ (Chandler et al. 1996; Gomez et al. 1997). According to Motte et al. (2001), TMC1A has an extended envelope and is recognized as a ‘true protostar’ not a ‘transition object’ nor an ‘edge-on view of T Tauri star’. Thus, we think that TMC1A is an archetype of Class I objects and an appropriate example to test our new method to assess the inclination. Observational data of TMC1A are taken from Myers et al. (1987), Kenyon et al. (1993), and Chandler et al. (1998).

Here, we adopt the new method to derive the inclination, and compare the estimate with the results by the full SED fitting. First, we attempt to estimate the inclination of the object using our new method. We evaluate first the observed luminosity $`L_{\mathrm{SED}}`$ by integrating the SED and obtain $`L_{\mathrm{SED}}=2.8L_{}`$, which is consistent with the results obtained by Myers et al. (1987) ($`L=2.4L_{}`$) and Kenyon & Hartmann (1995) ($`L=2.2L_{}`$). We next calculate $`f_\mathrm{L}`$ from $`L_{\mathrm{SED}}`$ and $`(\nu L_\nu )_{\mathrm{max}}`$, and obtain $`f_\mathrm{L}=2.52`$. From Fig. 4(a), this value indicates that inclination angle is about $`20^{}`$. This inclination is estimated with our new method only using the $`L_{\mathrm{SED}}`$ and the peak flux, without fitting all the physical parameters of the object.

Next, we try to obtain all the physical parameters of the object by fitting the SED and estimate the inclination angle. In practice, we fit the SED going through the following steps:

1. Set the outflow half opening angle $`\theta _{\mathrm{bp}}=20^{}`$.
2. Infer the luminosity of the central star $`L_{}`$ from the peak flux of the SED.
3. Infer the circumstellar mass $`M_{\mathrm{env}}+M_{\mathrm{disk}}`$ ($`M_{\mathrm{env}}\rho _1`$, $`M_{\mathrm{disk}}\mathrm{\Sigma }_1`$) from the flux in a frequency range lower than the frequency of the peak flux $`\nu _{\mathrm{peak}}`$.
4. Infer the ratio of the envelope mass to the disk mass and inclination from the shape of the SED in a frequency range higher than $`\nu _{\mathrm{peak}}`$.

The best fitted SED is displayed in Fig. 6 and adopted parameters are listed in Table 2. It is seen that our model can reproduce the observational data very well. In this procedure, the inclination of the object is estimated to be $`22^{}`$.

It is clear that the estimated values obtained by above two different methods agree quite well with each other. Since the estimation based on the full SED fitting is considered to be more reliable, it seems that the new spectrophotometric method provides an effective tool to assess the inclination.

However, it should be noted that the inclination angle obtained above does not agree with the values obtained by Kenyon et al. (1993), which is $`i60^{}`$, and Chandler et al. (1996), $`i=40^{}70^{}`$. Such discrepancy could be attributed to the method used to infer the inclination. Kenyon et al. (1993) used the spherically symmetric temperature distribution when they estimated the inclination with the SED. But the spherical temperature distribution can lead a significant difference in the resultant SED as shown in section 3.3. Chandler et al. (1996) estimated the inclination angle based on the shape of outflow lobe on the sky. They first restricted the inclination angle to be less than $`70^{}`$. Then, they placed the lower limit of the inclination to be $`40^{}`$, though the reason why they adopted the value was not given clearly.

### 4.3 Effects of Other Parameters

Here, we estimate uncertainties of the present method by changing parameters of the model. As shown in Fig. 4(a), an uncertainty of $`f_\mathrm{L}`$ is large if the half opening angle $`\theta _{\mathrm{bp}}`$ is thoroughly unknown. But, it is possible to pose a constraint for $`\theta _{\mathrm{bp}}`$ in the Class I phase by the observation of outflow lobes (Cabrit & Bertout 1986). For example, it is estimated that $`\theta _{\mathrm{bp}}=13^{}`$ \- $`19^{}`$ for Class I object TMC1 (IRAS04381+2540) and $`15^{}`$ \- $`21^{}`$ for TMC1A (IRAS04365+2535), respectively (Chandler et al. 1996). There are also some other results for estimation of $`\theta _{\mathrm{bp}}`$: $`\theta _{\mathrm{bp}}=22.5^{}`$ for L1448 (Bachiller et al. 1995) and B335 (Hirano et al. 1988; Chandler & Sargent 1993), and $`35^{}`$ for L1551 IRS5 (André et al. 1990). From these estimations, it seems possible to constrain $`\theta _{\mathrm{bp}}`$ to be $`10^{}`$ \- $`40^{}`$. This constraint would be justified by the fact that if $`\theta _{\mathrm{bp}}`$ becomes larger than this range, resultant SED exhibits double peaked feature in far-IR and near-IR to optical, but these SEDs are not classified Class I category. From Fig. 4(a), a typical error for determination of $`i`$ is $`\pm 10^{}`$ if we restrict ourselves to $`\theta _{\mathrm{bp}}=10^{}`$ \- $`40^{}`$.

Next, we examine to what extent the luminosity of the central star $`L_{}`$ affects the ratio $`f_\mathrm{L}`$. Fig. 7 shows calculated $`f_\mathrm{L}`$ as functions of the inclination $`i`$ with three different luminosities, $`0.1L_{}`$, $`1L_{}`$, and $`10L_{}`$. It is seen that when the luminosity of the central star is low ($`L_{}=0.1L_{}`$), $`f_\mathrm{L}`$ for the small inclinations decreases compared to the standard case ($`L_{}=1L_{}`$), but $`f_\mathrm{L}`$ for the large inclinations does not change. In contrast, when the luminosity is high ($`L_{}=10L_{}`$), $`f_\mathrm{L}`$ for the large inclinations becomes higher than the standard case, while $`f_\mathrm{L}`$ for the low inclinations remains the standard value. It is true that $`f_\mathrm{L}`$ is affected by the central star luminosity to some extent, but the change of the value in a range sensitive to the inclination is not large. Thus, when we estimate the inclination of an object using this new method and $`f_\mathrm{L}`$, the effect of the central star luminosity seems to be negligible.

Also, $`\rho _1`$ could affect the results. Fig. 8 shows the condition of $`f_\mathrm{L}=\mathrm{constant}`$ for each $`\rho _1`$ in the $`i\theta _{\mathrm{bp}}`$ plane. The flux at $`\nu >\nu _{\mathrm{peak}}`$ decreases with increase of $`\rho _1`$, because $`\tau _\mathrm{c}`$ is proportional to $`\rho _1`$ and the flux in this frequency range is proportional to $`\mathrm{exp}(\tau _\mathrm{c})`$. In contrast, the flux at $`\nu <\nu _{\mathrm{peak}}`$ increases with increase of $`\rho _1`$ because the total amount of dust grains is increased. Thus, $`f_\mathrm{L}`$ tends to decrease with increase of $`\rho _1`$, so that the line of $`f_\mathrm{L}=\text{constant}`$ is shifted above with increasing $`\rho _1`$ in the $`i\theta _{\mathrm{bp}}`$ plane. For $`\theta _{\mathrm{bp}}=10^{}40^{}`$, the inclination $`i`$ changes by $`10^{}`$ when $`\rho _1`$ changes by an order.

By contrast, $`\mathrm{\Sigma }_1`$ does not influence $`f_\mathrm{L}`$. Fig. 9 displays calculated $`f_\mathrm{L}`$ with different $`\mathrm{\Sigma }_1`$. It is easily seen that $`\mathrm{\Sigma }_1`$ does not change $`f_\mathrm{L}`$ significantly. This is because the circumstellar disk contribute little to the column density along the line of sight, on which $`f_\mathrm{L}`$ strongly depends.

Therefore, we conclude that the primary uncertainty of our method comes from the information of $`\theta _{\mathrm{bp}}`$. If $`\theta _{\mathrm{bp}}`$ can be estimated by any other means, then the error of inclination obtained by our method becomes roughly $`\pm 5^{}`$. If $`\theta _{\mathrm{bp}}`$ is unknown, then the error may be roughly $`\pm 15^{}`$.

## 5 Conclusions

Using 2D radiation transfer calculations, we have obtained the radiation fields for 2D axisymmetric protostar model, consistent with the central star, the circumstellar disk, and the envelope, and derived the SEDs of Class I objects. We have found that the ratio $`f_\mathrm{L}`$ between the emergent luminosity, $`L_{\mathrm{SED}}`$, and the peak flux in the SED, $`(\nu L_\nu )_{\mathrm{max}}`$, is a good indicator of the inclination angle of the object, because $`L_{\mathrm{SED}}`$ is sensitive to the inclination, while the peak flux is insensitive. Through the test with real data and the analysis on the other physical parameters, it has been shown that $`f_\mathrm{L}`$ is a robust tool to assess the inclination of a Class I object. Hence, $`f_\mathrm{L}`$ can provide a new spectrophotometric method to estimate the inclination angle of a Class I object. It is beneficial that, in this new method, only the luminosity $`L_{\mathrm{SED}}`$ and the peak flux $`(\nu L_\nu )_{\mathrm{max}}`$ are needed to estimate the inclination. The typical error of the method is roughly $`\pm 5^{}`$ if the half opening angle $`\theta _{\mathrm{bp}}`$ is known and roughly $`\pm 15^{}`$ if it is unknown.

The present method is applicable for a great deal of data provided by future missions of infrared, sub-millimeter, and millimeter observatories, e.g., SIRTF, ASTRO-F, and ALMA. These missions are expected to reveal numerous protostar candidates. Then, it may be possible to carry out precise statistical study of protostars. The present spectrophotometric method could be a powerful tool in the statistical study.

We are grateful to Dr. N. Kikuchi for providing us the original numerical code. Numerical calculations were carried out with the facilities in Center for Computational Physics, University of Tsukuba. This work is supported in part by the Grant-in-Aid of the JSPS (TN 12740118 and 10147105; MU 11640255).


format: text

text:
# Identical Particle Interferometry at STAR

## 1 Introduction

Two-particle intensity interferometry (HBT) is a useful tool to study the space-time geometry of the particle-emitting source in heavy ion collisions . It also contains dynamical information that can be explored by studying the transverse momentum ($`p_T`$) dependence of the apparent source size . In non-central collisions information on the anisotropic shape of the pion-emitting region can be extracted by measuring two-pion correlation functions as a function of emission angle with respect the reaction plane .

In this paper we present two-pion correlation systematics as a function of the collision energy $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}`$, transverse mass ($`m_T=\sqrt{p_T^2+m^2}`$), multiplicity, and emission angle with respect to the reaction plane in Au+Au collisions produced by the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory. We also present the first significant HBT measurement for neutral kaons from heavy ion collisions; besides carrying important information about the dynamics of strange particles, these measurements provide valuable cross-checks on charged kaon HBT, as complication due to various corrections (e.g. two-track efficiency and Coulomb correction) are quite different for neutral particles.

Experimentally, two-particle correlations are studied by constructing the correlation function $`C_2`$(q) = A(q)/B(q). Here A(q) is the measured distribution of the momentum difference q = $`\text{p}_1`$ \- $`\text{p}_2`$ for pairs of particles from the same event, and B(q) is the corresponding distribution for pairs of particles from different events.

## 2 Experimental Details

For this analysis we selected events with a collision vertex position within $`\pm `$30 cm measured from the center of the 4 m long STAR Time Projection Chamber (TPC), and we mixed events only if their longitudinal primary vertex positions were no farther apart than 6 cm. We divided our sample into three centrality bins, where the centrality was characterized according to the measured multiplicity of charged particles at midrapidity. The three centrality bins correspond to 0-10$`\%`$ (central), 10-30$`\%`$ (midcentral), and 30-70$`\%`$ (peripheral) of the total hadronic cross section. Charged pions were identified by correlating their specific ionization in the gas of the TPC with their measured momentum. Neutral kaons were identified via topological methods .

For charged pions, the effects of track-splitting (reconstruction of a single track as two tracks) and track-merging (two tracks with similar momenta reconstructed as a single track) were eliminated as described in . We applied to each background pair a Coulomb correction  corresponding to a spherical Gaussian source of 5 fm radius.

In the neutral particle analysis, the possibility of a single $`K_s^0`$ being correlated with itself was eliminated by requiring that a pair of $`K_s^0`$’s have unique daughters, and that their decay positions were spatially well-separated.

## 3 Pion HBT versus $`p_T`$ and centrality

The three-dimensional correlation functions were generated. The relative momentum q was measured in the longitudinal co-moving system (LCMS) frame, and decomposed according to the Pratt-Bertsch  ”out-side-long” parametrization. The correlation functions were fit with: C($`q_o`$$`q_s`$$`q_l`$) = 1 + $`\lambda `$ $`\mathrm{exp}`$(-$`R_o^2`$ $`q_o^2`$ \- $`R_s^2`$ $`q_s^2`$ \- $`R_l^2`$ $`q_l^2`$).

The effect of the single-particle momentum resolution ($`\delta `$p/p $``$ 1$`\%`$ for the particles under study) induces systematic underestimation of the HBT parameters. Using an iterative procedure , we corrected our correlation functions for finite resolution effects. The correction due to the uncertainty on the removal of the artificial reduction of the HBT parameters associated with the anti-merging cut is still being finalized; we estimate systematic errors on HBT radii at about 1 fm.

Fig.2 shows the $`m_T`$ dependence of the source parameters for negative pions at three centrality bins. The three radii increase with increasing centrality and $`R_l`$ varies similarly to $`R_o`$, $`R_s`$; for $`R_o`$ and $`R_s`$ this increase may be attributed to the geometrical overlap of the two nuclei. The extracted radii rapidly decrease as a function of $`m_T`$, which is an indication of transverse flow ; interestingly the $`m_T`$ dependence is independent of centrality.

While slightly larger at 200 GeV, the transverse homogeneity lengths are similar for both energies (Figs.2 (a) and (b)). These results suggest a steeper fall-off with $`m_T`$, which might indicate increased flow; finalized corrections are necessary before drawing such conclusions. $`R_o`$ falls somewhat steeper than $`R_s`$ with $`m_T`$ and $`R_o/R_s1`$ (Fig. 2(d)) which indicates short emission duration, $`\mathrm{\Delta }`$t $``$ 1-2 fm/c, in a blast wave fit .

The longitudinal radius falls along the same curve at both energies (Fig. 2(c)). Assuming boost-invariant longitudinal flow, we can extract an evolution timescale, by using a simple fit  (solid lines Fig. 2(c)): $`R_l=<t_{fo}>\sqrt{\frac{T_k}{m_T}}`$. Taking $`T_k`$ = 110 MeV we got $`<t_{fo}>`$ $``$ 10 fm/c for central events and $`<t_{fo}>`$ $``$ 7.6 fm/c for peripheral events. Hence, the evolution time, in addition to the emission duration, is quite short.

## 4 Pion HBT with respect the reaction plane

Azimuthally-sensitive HBT, correlated event-by-event with the reaction plane, provides information about the coordinate-space anisotropies of the emitting source .

We calculated the three-dimensional correlation functions with fixed pair angle with respect to the $`2^{nd}`$-order event plane (determined from elliptic flow) $`\varphi `$ = $`\mathrm{}(\text{K}_{},\text{b})`$, where $`\text{K}_{}`$ = $`(\text{p}_1+\text{p}_2)_{}`$ is the total momentum of the pair perpendicular to the beam. The correlation functions for each $`\varphi `$ bin were fit with the standard Gaussian parametrization C(q,$`\varphi `$) = 1 + $`\lambda (\varphi )\mathrm{exp}[q_iq_jR_{ij}^2(\varphi )]`$, in this case the cross-term $`R_{os}^2`$ is also relevant.

Fig.2 shows the HBT parameters plotted as a function of $`\varphi `$ at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}`$=130 GeV. Data were corrected for event plane resolution and merging systematics. $`R_o`$ and $`R_s`$ show significant equal and opposite $`2^{nd}`$-order oscillations with phases indicating an out of plane extended source geometry. This would suggest a fast evolution, since positive $`v_2`$ will push the source from an out-of-plane extended configuration (present in the entrance channel and caused by the partial overlap of the colliding nuclei) to in-plane extended. For comparison, the lines represent a fit with the ”blast wave model” , corresponding to a transversely deformed (5$`\%`$ extended out-of-plane) source geometry with an oscillating flow field, and a short (2 fm/c) emission duration. This model successfully reproduces the details of elliptic flow measured by STAR .

Preliminary results for $`\pi ^{}`$ and $`\pi ^+`$ from peripheral Au+Au collisions at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}`$=200 GeV confirm the oscillations measured at 130 GeV and the higher statistics of these measurements will allow $`m_T`$ and centrality systematic analysis.

## 5 $`K_s^0`$ Interferometry

The good efficiency of the TPC for finding and reconstructing short-lived neutral kaons gives us the opportunity to study HBT correlations for $`K_s^0`$. This analysis should allow us to extend our systematics to higher $`p_T`$. Using the topology of the decay $`K_s^0\pi ^+\pi ^{}`$ we reconstructed $``$3.8 $`K_s^0`$ candidates per event with $`<m_T>`$ 1.12 GeV/c2. Fig.3 shows the one-dimensional correlation function which exhibits a promising low-Q correlation signal. Fitting with the functional form C(Q) = 1 + $`\lambda \mathrm{exp}(Q^2R^2)`$ yields $`\lambda =0.76\pm 0.29`$ and $`R=5.75\pm 1`$ fm. If these results persist as the analysis is finalized, it is interesting to speculate on possible causes of such a large homogeneity length at this high $`m_T`$.

## 6 Conclusion

We have presented identical meson interferometry results for Au+Au collisions at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}`$=200 GeV. With respect to multiplicity and $`m_T`$ dependencies, pion HBT radii are very similar to results reported at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}`$=130 GeV. Our results indicate that both the evolution timescale (as measured by the $`m_T`$ dependence of $`R_l`$ and the out-of-plane extended source indicated by the azimuthally-sensitive HBT) and the emission duration (probed by comparing $`R_o`$ to $`R_s`$) are surprisingly fast. The large acceptance and excellent tracking of the TPC has allowed for the first solid measurement of $`K_s^0`$-$`K_s^0`$ correlations. While still preliminary, these suggest a surprisingly large homogeneity length at high $`m_T`$.


text:
# Near Horizon Geometry of Extreme Black Holes and Colliding Waves

## Abstract

The correspondence between black holes and colliding waves extends to cover the near horizon geometry of rotating black holes and colliding waves with cross polarization. Extreme Kerr and Kerr-Newman geometries are given as examples.

It is believed that matter collisions at speeds almost the speed of light can produce mini black holes even in a laboratory. This possibility directs attentions naturally on collisions of waves in general relativity. Local equivalence (isometry) between two different types of spacetimes, such as black holes (BHs) and colliding plane waves (CPWs) was pointed out first by Chandrasekhar and Xanthopoulos CX . This isometry covered only a portion of the BH geometry, namely region between horizons and the interaction region of CPWs. This concerns only those CPWs which create horizon instead of an all encompassing space-like singularity. Such examples of CPWs are very rare so far. The prototype is the solution of CX which describes collision of impulse gravitational waves accompanied by shock gravitational waves. As a second example we quote colliding electromagnetic (em) shock waves of Bell and Szekeres (BS) , which is transformable to the Bertotti-Robinson (BR) solution . This latter is important from its correspondence to conformal field theory and the underlying anti-de Sitter (AdS) structure with a group $`SL(2,R)\times SO(3).`$ From geometrical point of view BR represents the throat region of an extreme Reissner-Nordstrom BH and its equivalence to colliding em waves is well-known .

In this letter we show that the connection between BHs and CPWs  is not restricted only to the local isometry of CX, but seems to be more deep rooted. This we do by showing that the near horizon geometry of the Extreme Kerr and Kerr-Newman BHs are identical with particular CPWs.

The Kerr metric in the Boyer-Linquist coordinates is given by

$$ds^2=e^{2\nu }d\stackrel{~}{t}^2e^{2\psi }\left(d\stackrel{~}{\varphi }\omega d\stackrel{~}{t}\right)^2\rho ^2\left(\frac{d\stackrel{~}{r}}{\stackrel{~}{\mathrm{\Delta }}}+d\theta ^2\right)$$
(1)
where

$`e^{2\nu }`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}\rho ^2}{\left(\stackrel{~}{r}^2+a^2\right)^2\stackrel{~}{\mathrm{\Delta }}a^2\mathrm{sin}^2\theta }}`$ (2)
$`e^{2\left(\psi +\nu \right)}`$ $`=`$ $`\stackrel{~}{\mathrm{\Delta }}\mathrm{sin}^2\theta `$ (3)
$`\omega e^{2\nu }`$ $`=`$ $`{\displaystyle \frac{2Ma\stackrel{~}{r}}{\stackrel{~}{\mathrm{\Delta }}\rho ^2}}`$ (4)
$`\rho ^2`$ $`=`$ $`\stackrel{~}{r}^2+a^2\mathrm{cos}^2\theta `$ (5)
$`\stackrel{~}{\mathrm{\Delta }}`$ $`=`$ $`\left(\stackrel{~}{r}r_+\right)\left(\stackrel{~}{r}r_{}\right)`$ (6)
$`r_+`$ $`=`$ $`M+\sqrt{M^2a^2}`$ (7)
$`r_{}`$ $`=`$ $`M\sqrt{M^2a^2}`$ (8)
In the extreme case we have $`M=a,`$ so that angular momentum is $`J=M^2,`$ the horizon area becomes $`A=8\pi M^2`$ while $`\stackrel{~}{\mathrm{\Delta }}=\left(\stackrel{~}{r}M\right)^2.`$ The near horizon geometry is described by setting

$`\stackrel{~}{r}`$ $`=`$ $`M+\lambda r`$ (9)
$`\stackrel{~}{t}`$ $`=`$ $`{\displaystyle \frac{t}{\lambda }}`$ (10)
$`\stackrel{~}{\varphi }`$ $`=`$ $`\varphi +{\displaystyle \frac{t}{2M\lambda }}`$ (11)
and taking the limit $`\lambda 0.`$ One obtains (after shifting also $`r\frac{1}{r}`$ ) the geometry

$$ds^2=\frac{\left(1+\mathrm{cos}^2\theta \right)}{2}\left[\left(\frac{dt}{r_0r}\right)^2\left(\frac{r_0dr}{r}\right)^2\left(r_0d\theta \right)^2\right]\frac{2\left(r_0\mathrm{sin}\theta \right)^2}{\left(1+\mathrm{cos}^2\theta \right)}\left(d\varphi \frac{dt}{rr_0^2}\right)^2$$
(12)
in which $`r_0^2=2M^2`$ . We set further, for simplicity $`r_0^2=2M^2=1.`$ Although BR metric represents the throat region of an extreme Reissner-Nordstrom BH the metric (4) due to itscross term signals more of a BH property. It can easily be checked that the Killing vector $`_t`$ becomes spacelike for $`\mathrm{cos}2\theta <4\sqrt{3}7.`$ We show now that this near horizon geometry of an extreme Kerr BH coincides with the colliding wave metric of CX . The CX metric describing collision of impulsive gravitational waves accompanied by shock waves is given by

$$ds^2=X\left(\frac{d\eta ^2}{\mathrm{\Delta }}\frac{d\mu ^2}{\delta }\right)\mathrm{\Delta }\delta \frac{X}{Y}dy^2\frac{Y}{X}\left(dxq_2dy\right)^2$$
(13)
where

$`X`$ $`=`$ $`\left(1p\eta \right)^2+q^2\mu ^2`$ (14)
$`Y`$ $`=`$ $`p^2\mathrm{\Delta }+q^2\delta `$ (15)
$`q_2`$ $`=`$ $`{\displaystyle \frac{2q\left(1\eta \right)\left(p\mu ^2+p\eta +\mu ^21\right)}{\left(1+p\right)\left(1p^2\eta ^2q^2\mu ^2\right)}}`$ (16)
$`\mathrm{\Delta }`$ $`=`$ $`1\eta ^2`$ (17)
$`\delta `$ $`=`$ $`1\mu ^2`$ (18)
in which the constants of second polarization $`(p,q)`$ are constrained by

$$p^2+q^2=1$$
(19)
We note that the metric function $`q_2`$ is determined up to an additive constant. This metric admits a horizon at $`\eta =1`$ and is free of curvature singularities. Complete analytic extension of this space time was given in detail by CX . In particular the extended spacetime has timelike singularities on hyperbolic arcs reminiscent of the ring singularity in the Kerr space. Now let us choose $`p=0,q=1`$ so that the metric (5) becomes

$$ds^2=X\left(\frac{d\eta ^2}{\mathrm{\Delta }}\frac{d\mu ^2}{\delta }\mathrm{\Delta }dy^2\right)\frac{4Y}{X}\left(d\stackrel{~}{x}\eta dy\right)^2$$
(20)
where

$`X`$ $`=`$ $`1+\mu ^2`$ (21)
$`Y`$ $`=`$ $`1\mu ^2`$ (22)
and we have introduced the new coordinate

$`d\stackrel{~}{x}={\displaystyle \frac{1}{2}}dx+dy`$

Next we apply the transformation

$`\eta `$ $`=`$ $`{\displaystyle \frac{1}{2r}}\left(r^2t^2+1\right)`$ (23)
$`\mathrm{tanh}y`$ $`=`$ $`{\displaystyle \frac{1}{2t}}\left(r^2t^21\right)`$ (24)
$`\mu `$ $`=`$ $`\mathrm{cos}\theta `$ (25)
$`\stackrel{~}{x}`$ $`=`$ $`\varphi {\displaystyle \frac{1}{2}}\mathrm{ln}\left[{\displaystyle \frac{\left(r+t\right)^21}{\left(rt\right)^21}}\right]`$ (26)
which results in

$$ds^2=\left(1+\mathrm{cos}^2\theta \right)\left(\frac{dt^2dr^2}{r^2}d\theta ^2\right)\frac{4\mathrm{sin}^2\theta }{1+\mathrm{cos}^2\theta }\left(d\varphi \frac{1}{r}dt\right)^2$$
(27)
This metric is (up to a scale factor of $`\frac{1}{2}`$ ) the line element (4) for the geometry of near horizon limit of an extreme Kerr metric. We note that by scaling an overall constant $`Q^2`$ can be introduced in front of this metric. This amounts, in the transformation (10) to taking $`\stackrel{~}{x}Q\stackrel{~}{x}`$ and $`yQy`$. We add also that the choice $`p=0,q=1`$ pushes the timelike singularities on hyperbolic arcs to infinity .

The Kerr-Newman metric has the same form as (1) except that $`\stackrel{~}{\mathrm{\Delta }}=\stackrel{~}{r}^22M\stackrel{~}{r}+a^2+q^2`$, where $`q`$ is the electric charge and the $`2M\stackrel{~}{r}`$ factor in $`\omega `$ is replaced by $`\stackrel{~}{r}^2+a^2\stackrel{~}{\mathrm{\Delta }}`$ . The extreme Kerr-Newman case corresponds to $`M^2=a^2+q^2`$ and the horizon is at $`\stackrel{~}{r}=M`$ with area $`4\pi \left(M^2+a^2\right)`$ . Expanding $`\omega `$ to first order in $`\stackrel{~}{r}M`$ one obtains

$`\omega `$ $`=`$ $`{\displaystyle \frac{a}{r_0^2}}{\displaystyle \frac{2Ma}{r_0^4}}\left(\stackrel{~}{r}M\right)`$ (28)
$`r_0^2`$ $`=`$ $`M^2+a^2`$ (29)
the near horizon limit takes the form $`\left(\text{with }r\frac{1}{r}\right)`$

$`ds^2`$ $`=`$ $`\left(1{\displaystyle \frac{a^2\mathrm{sin}^2\theta }{r_0^2}}\right)\left[\left({\displaystyle \frac{dt}{r_0r}}\right)^2\left({\displaystyle \frac{r_0dr}{r}}\right)^2\left(r_0d\theta \right)^2\right]`$ (31)
$`r_0^2\mathrm{sin}^2\theta \left(1{\displaystyle \frac{a^2\mathrm{sin}^2\theta }{r_0^2}}\right)^1\left(d\varphi +{\displaystyle \frac{2Ma}{rr_0^4}}dt\right)^2`$
We proceed now to show that this corresponds to a metric of colliding cross polarized em shock waves . In the notation of CX this metric is given by

$$ds^2=F\left(\frac{d\eta ^2}{\mathrm{\Delta }}\frac{d\mu ^2}{\delta }\right)\mathrm{\Delta }\delta dy^2\frac{\delta }{F}\left(dxq_0\eta dy\right)^2$$
(32)
where

$`2F`$ $`=`$ $`1+p_0+\left(p_01\right)\mu ^2`$ (33)
$`p_0`$ $`=`$ $`\sqrt{1+q_0^2}=const.`$ (34)
In terms of the null coordinates $`(u,v)`$ we have $`\eta =sin(au+bv)`$ and $`\mu =sin(aubv)`$, where $`(a,b)`$ are the constant em parameters. For $`q_0=0`$ $`(p_0=1)`$ this solution reduces to the BS solution of colliding em waves with linear polarization. Applying a similar transformation as (10) (with $`\stackrel{~}{x}=x`$ now) casts our metric into

$$ds^2=F(\theta )\left(\frac{dt^2dr^2}{r^2}d\theta ^2\right)\frac{\mathrm{sin}^2\theta }{F(\theta )}\left(d\varphi \frac{dt}{r}\right)^2$$
(35)
where

$$2F\left(\theta \right)=1+p_0+\left(p_01\right)\mathrm{cos}^2\theta $$
(36)
Comparing the metrics (13) and (16) we observe that upon choosing $`q_0=Ma`$ with $`M^2=a^2+2`$ (which corresponds to a charge of $`\sqrt{2}`$ ) and scalings $`\varphi r_0^2\varphi `$, $`t\frac{t}{r_0^2}`$ they become identical. This proves that the near horizon geometry of an extreme Kerr-Newman BH is obtained as a result of two colliding em waves with cross polarization.

A general proof, namely, that the extreme geometry of any BH is represented by a CPW remains open. Generic representations of both classes ( of BHs and CPWs ) may provide a general proof. We are satisfied here only with the two simplest but most important cases. A general proof on the other hand is reminiscent of a theorem by Penrose : Any spacetime admits a plane wave limit. Here it is not any spacetime but it is only the near horizon geometry of an extreme BH. And it is not a plane wave limit but it is a colliding wave limit

###### Acknowledgement 1

I wish to thank Ozay Gurtug for helpful discussions.


text:
# 1 Parameter area which we can explore through the asymmetry 𝐴_{ℓℓ}. We can confirm this asymmetry to be non-zero at 1𝜎, 2𝜎 and 3𝜎 level when the parameters "Re"(𝐷_{𝛾,𝑍}) giving 𝐶𝑃 violation in the 𝑡𝑡̄𝛾/𝑍 couplings and "Re"(𝑓^𝑅₂-𝑓̄^𝐿₂) giving 𝐶𝑃 violation in the 𝑡𝑏𝑊 coupling are outside the two solid lines, dashed lines and dotted lines respectively. Unfortunately there is some area where two contributions from the production and decay vertices cancel each other and we get little information.

TOKUSHIMA Report (hep-ph/0210224)

A New Decoupling Theorem in Top-Quark Physics ***Talk at International Workshop on Physics and Experiments with Future Electron-Positron Linear Colliders (LCWS2002), August 26 - 30, 2002, Jeju-island, Korea.
This work is based on collaboration with B. Grzadkowski.

Zenrō HIOKI$`^{}`$\*)\*)\*)E-mail address: hioki@ias.tokushima-u.ac.jp

Institute of Theoretical Physics, University of Tokushima

Tokushima 770-8502, JAPAN

ABSTRACT

Angular distribution of a secondary particle in top-quark production/decay is studied in a simple and general manner. It is shown that the distribution does not depend on any possible anomalous top-quark-decay interactions whatever the production mechanism is when certain well-justified conditions are satisfied. Some analyses using the final-state lepton are presented as an example of its application.

Top quark possesses a huge mass, which is very close to the electroweak (EW) scale. Therefore studying its property in various aspects will give us valuable information on the EW symmetry breakdown and consequently the origin of particle masses. Since future high-energy accelerators like NLC/LHC are expected to work as top-quark factories, a lot of attention has been paid to explore their productions (for a review, see Ref. and the reference list there).

We could classify those analyses into two categories: model-dependent one and model-independent one. In the former approach, we can perform precise calculations and compare the results with corresponding experimental data at high precision. However if the nature does not choose the model, we have to look for another candidate. On the other hand, the latter offers us a “no-lose game” in a sense since we could always get some information anyway, but we have to introduce a lot of free parameters (form factors) to describe the interactions as generally as possible. Thus it is not that easy to achieve a high precision. Therefore these two approaches have both good and bad points respectively, and should work complementary to each other.

During investigating $`e^+e^{}t\overline{t}\mathrm{}^\pm \mathrm{}`$ model-independently, we found that the angular distribution of the leptons $`\mathrm{}^\pm `$ in open-top region is not sensitive to modification of the standard $`tbW`$ decay vertex . The same conclusion was also reached by Rindani  through an independent calculation. Because of the above-mentioned reason, a distribution insensitive to a certain class of nonstandard form factors is obviously a big advantage in model-independent analyses.

This discovery was rather accidental and unexpected, but after that we have succeeded to show in Ref. that this phenomenon holds quite generally and what we found first was one typical example of it. Therefore we now believe it deserves to be called a new “Decoupling Theorem”, on which I would like to report here.

Let us consider a general top-quark-production process $`1+2t+\mathrm{}`$ followed by a decay $`tf+\mathrm{}`$, where $`f`$ denotes the secondary particle that we are interested in. Since the ratio of the top-quark width $`\mathrm{\Gamma }_t`$ to its mass $`m_t`$ is of the order of $`10^2`$, we can safely adopt the narrow-width approximation for the decaying top and apply the Kawasaki-Shirafuji-Tsai formula  in order to determine the $`f`$ distribution except in the threshold region♯1♯1♯1In the threshold region, this formula might be no longer valid due to large corrections. For example, non-factorizable QCD corrections appear at the level of 10 % in $`e\overline{e}t\overline{t}\mathrm{}^\pm X`$.:

$$\frac{d\sigma }{d\stackrel{~}{𝒑}_f}=2B_f𝑑\stackrel{~}{𝒑}_t\frac{d\sigma }{d\stackrel{~}{𝒑}_t}(s_t=n)\frac{1}{\mathrm{\Gamma }}\frac{d\mathrm{\Gamma }}{d\stackrel{~}{𝒑}_f}.$$
(1)
Here $`d\stackrel{~}{𝒑}`$ denotes the Lorentz-invariant phase-space element $`d𝒑/[(2\pi )^32p^0]`$, $`d\mathrm{\Gamma }/d\stackrel{~}{𝒑}_f`$ is the spin-averaged top-quark width to $`f+\mathrm{}`$, $`B_f\mathrm{\Gamma }/\mathrm{\Gamma }_t`$, and $`d\sigma (s_t=n)/d\stackrel{~}{𝒑}_t`$ is the single-top-quark inclusive cross section with the top-spin vector $`s_t`$ being replaced with the so-called “effective polarization vector” $`n`$, which is a function of $`p_{t,\mathrm{}}`$ and carries important information on the top decay. Its definition is a bit complicated, but it can be expressed as

$$n=\alpha ^f\left(\frac{m_t}{p_tp_f}p_f\frac{p_t}{m_t}\right)$$
(2)
thanks to Lorentz covariance, where $`\alpha ^f`$ is a real parameter constrained as $`|\alpha ^f|1`$.

The angular distribution of $`f`$ is obtained by integrating eq.(1) over $`E_f`$:

$$\frac{d\sigma }{d\mathrm{\Omega }_f}=\frac{B_f}{(2\pi )^3}𝑑E_fE_f𝑑\stackrel{~}{𝒑}_t\frac{d\sigma }{d\stackrel{~}{𝒑}_t}(s_t=n)\frac{1}{\mathrm{\Gamma }}\frac{d\mathrm{\Gamma }}{d\stackrel{~}{𝒑}_f}.$$
(3)
Here the polarization-vector $`n`$ may in general depend on $`E_f`$ and the integration over $`E_f`$ on the right-hand side cannot be performed any further without knowing explicit form of $`d\sigma (s_t=n)/d\stackrel{~}{𝒑}_t`$. However, if the vector $`n`$ is free from $`E_f`$, we can perform the $`E_f`$ integration independently of the production mechanism since $`d\sigma /d\stackrel{~}{𝒑}_t`$ can depend on $`E_f`$ only through $`n`$:

$$\frac{d\sigma }{d\mathrm{\Omega }_f}=2B_f𝑑\stackrel{~}{𝒑}_t\frac{d\sigma }{d\stackrel{~}{𝒑}_t}(s_t=n)\frac{1}{\mathrm{\Gamma }}\frac{d\mathrm{\Gamma }}{d\mathrm{\Omega }_f}.$$
(4)
Since $`d\mathrm{\Gamma }/d\stackrel{~}{𝒑}_f`$ in eq.(1) is the unpolarized top width as was explained, the resultant angular distribution is isotropic in the top-quark rest frame. Therefore its form in the Lab frame is fully determined by the Lorentz transformation connecting these two frames as

$$\frac{d\mathrm{\Gamma }}{d\mathrm{cos}\theta _{tf}}=\frac{1\beta ^2}{(1\beta \mathrm{cos}\theta _{tf})^2}\frac{d\mathrm{\Gamma }^{}}{d\mathrm{cos}\theta ^{}},$$
(5)
where $`d\mathrm{\Gamma }^{}/\mathrm{d}\mathrm{cos}\theta ^{}`$ is the constant distribution defined in the top-quark rest frame and $`\theta _{tf}`$ means the angle between $`𝒑_t`$ and $`𝒑_f`$ in the Lab frame. Thus we obtain

$$\frac{d\sigma }{d\mathrm{\Omega }_f}=\frac{1}{2\pi }B_f𝑑\stackrel{~}{𝒑}_t\frac{d\sigma }{d\stackrel{~}{𝒑}_t}(s_t=n)\frac{1\beta ^2}{(1\beta \mathrm{cos}\theta _{tf})^2}.$$
(6)
Note that there are only two possible ways that the structure of the top-quark-decay vertices could influence the distribution:
i) through the width $`d\mathrm{\Gamma }/d\stackrel{~}{𝒑}_f`$, ii) through the effective polarization vector $`n`$.
Therefore we conclude that if the polarization vector $`n`$ depends neither on $`E_f`$ nor on anomalous top-quark-decay vertices, the angular distribution $`d\sigma /d\Omega _f`$ is not altered by those anomalous vertices except for possible trivial modification of the branching ratio $`B_f`$. Furthermore, if we focus on the single standard-decay channel $`tbWb\mathrm{}\nu _{\mathrm{}}`$, even that dependence disappears.

Now the question is: Is $`n`$ really $`E_f`$ independent? In the framework of the standard model $`tbWb\mathrm{}\nu _{\mathrm{}}`$ is practically the only decay mode of the top quark, and it was found in Ref. that

$$\alpha ^\mathrm{}^+=1\mathrm{and}\alpha ^b=(2M_W^2m_t^2)/(2M_W^2+m_t^2).$$
(7)
This result shows that the form of $`n`$ vector is process-dependent but does not depend on $`E_f`$ indeed at least within the SM if $`m_f`$ can be neglected.

This process will still be the main decay mode even after taking general couplings into account, unless we unrealistically change the top-interaction structure. Therefore, we assumed the following most general covariant $`tbW`$ coupling and calculated $`n`$ explicitly:

$$\mathrm{\Gamma }^\mu \overline{u}(p_b)\left[\gamma ^\mu (f_1^LP_L+f_1^RP_R)\frac{i\sigma ^{\mu \nu }k_\nu }{M_W}(f_2^LP_L+f_2^RP_R)\right]u(p_t),$$
(8)
where $`P_{L/R}=(1\gamma _5)/2`$ and $`k`$ is the momentum of $`W`$. As a result, we have observed that $`n(f=\mathrm{}^+)`$ remains unchanged while $`n(f=b)`$ receives corrections. In those calculations, all the fermions except $`t`$ and $`b`$ were treated as massless,♯2♯2♯2Some part of the correction to $`n(f=b)`$ does not vanish even for $`m_b=0`$. the narrow-width approximation was adopted also for the decaying $`W`$, and only the $`[`$SM$`]`$-$`[`$non-SM$`]`$ interference terms were taken into account. Thus, any anomalous $`tbW`$ interactions decouple from the leptonic angular distribution within this approximation.

Let me show some application of this “Decoupling Theorem”. We introduced in Ref. the following $`CP`$-violating asymmetries:

$$𝒜_{CP}(\theta )=\left[\frac{d\sigma ^+(\theta )}{d\mathrm{cos}\theta }\frac{d\sigma ^{}(\pi \theta )}{d\mathrm{cos}\theta }\right]/\left[\frac{d\sigma ^+(\theta )}{d\mathrm{cos}\theta }+\frac{d\sigma ^{}(\pi \theta )}{d\mathrm{cos}\theta }\right],$$
(9)
$$𝒜_{CP}=\frac{{\displaystyle _1^0}d\mathrm{cos}\theta {\displaystyle \frac{d\sigma ^+(\theta )}{d\mathrm{cos}\theta }}{\displaystyle _0^{+1}}d\mathrm{cos}\theta {\displaystyle \frac{d\sigma ^{}(\theta )}{d\mathrm{cos}\theta }}}{{\displaystyle _1^0}d\mathrm{cos}\theta {\displaystyle \frac{d\sigma ^+(\theta )}{d\mathrm{cos}\theta }}+{\displaystyle _0^{+1}}d\mathrm{cos}\theta {\displaystyle \frac{d\sigma ^{}(\theta )}{d\mathrm{cos}\theta }}}.$$
(10)
Of course they are a pure measure of the $`CP`$-violating anomalous top productions. These lepton asymmetries are quite in contrast to the following asymmetry

$$A_{\mathrm{}\mathrm{}}\frac{{\displaystyle _{x<\overline{x}}𝑑x𝑑\overline{x}\frac{d^2\sigma }{dxd\overline{x}}}{\displaystyle _{x>\overline{x}}𝑑x𝑑\overline{x}\frac{d^2\sigma }{dxd\overline{x}}}}{{\displaystyle _{x<\overline{x}}𝑑x𝑑\overline{x}\frac{d^2\sigma }{dxd\overline{x}}}+{\displaystyle _{x>\overline{x}}𝑑x𝑑\overline{x}\frac{d^2\sigma }{dxd\overline{x}}}}$$
(11)
introduced in Ref. using the $`\mathrm{}^\pm `$ energy correlation $`d^2\sigma /dxd\overline{x}`$, where $`x`$ and $`\overline{x}`$ are normalized energies of $`\mathrm{}^+`$ and $`\mathrm{}^{}`$ respectively. Generally this is also an asymmetry very sensitive to nonstandard $`CP`$ violation in the top-quark interactions. However, if we have no luck and two contributions from the production and decay vertices cancel each other, we end up with getting little information as found in Fig.1. This comparison lightens the outstanding feature of $`𝒜_{CP}(\theta )`$ and $`𝒜_{CP}`$ more clearly. This is however no longer a defect of $`A_{\mathrm{}\mathrm{}}`$. If we get to know the top-production mechanism through $`𝒜_{CP}/𝒜_{CP}(\theta )`$, then we can thereby explore the decay mechanism via $`A_{\mathrm{}\mathrm{}}`$. This would never be possible if we had $`𝒜_{CP}/𝒜_{CP}(\theta )`$ alone. That is, our decoupling theorem turns the defect of $`A_{\mathrm{}\mathrm{}}`$ to an advantage!

Finally let me give some comments: Our derivation of the angular distribution (6) applies to any top-quark production process, including pair and single productions at both $`e^+e^{}`$ and hadronic machines (or $`\gamma \gamma `$ collisions enabled by laser-electron/positron backward scatterings). For $`e^+e^{}`$ collisions the absolute value of top-quark momentum is fixed by $`\beta ^2=14m_t^2/s`$ and eq.(6) reduces to

$$\frac{d\sigma }{d\mathrm{\Omega }_f}=\frac{2m_t^2}{\pi s}B_f𝑑\mathrm{\Omega }_t\frac{d\sigma }{d\mathrm{\Omega }_t}(s_t=n)\frac{1}{(1\beta \mathrm{cos}\theta _{tf})^2}.$$
(12)
This agrees with the one derived by Arens and Sehgal within the SM , where they pointed out that the lepton angular distribution can be used as a good spin analyzer of the parent top quark. Our theorem assures that their conclusion is not affected by the anomalous $`tbW`$ couplings. On the other hand, the distribution in the CM frame of hadron-hadron collisions has some additional factors since the hadron CM frame and the parton CM frame are different from each other and they are connected through Lorentz transformation. However, any Lorentz boost can never produce anomalous-decay-parameter dependence. So, if $`d\sigma /d\mathrm{cos}\theta `$ in the parton-CM frame is free from the non-SM form factors, then the one in the hadron-CM frame is also free from them. Consequently, our decoupling theorem holds in hadron-hadron collisions, too.

In summary, we have investigated the angular distribution of a secondary particle $`f`$ in processes like $`1+2t+\mathrm{}`$ followed by $`tf+\mathrm{}`$ neglecting the $`f`$ mass and applying the narrow-width approximation for the decaying top. It has been clarified that if the effective polarization vector $`n`$ contains neither non-SM top-quark couplings nor $`E_f`$ the whole angular distribution of $`f`$ has no non-SM top-quark-decay contributions. We then showed that this is realized in one of the most significant cases $`f=\mathrm{}^+`$ within our approximation. We expect that NLC/LHC will be able to give certain clear statements about the top-quark interactions through measurements of the secondary lepton distributions.


text:
# Dileptons and Photons from Coarse–Grained Microscopic Dynamics and Hydrodynamics Compared to Experimental Data

## 1 Introduction

We focus on central Pb-Pb and Pb-Au collisions at 158 GeV per nucleon at the CERN SPS. A coarse–graining of 100 UrQMD simulation events was performed. This gives flow fields of temperature, baryon chemical potential, and fluid velocity that can be compared to the flow fields of relativistic hydrodynamics. The initial conditions and final freezeout temperature in the hydrodynamic calculations were previously adjusted to represent the rapidity and transverse momentum distributions of hadrons so there are no new parameters introduced here. The coarse–grained version of UrQMD does require a freezeout temperature; a value of 125$`\pm `$10 MeV reproduces not only the rapidity distribution very well (shown) but also the transverse momentum distributions of negative hadrons and the net proton number (not shown). We then compare to the measured photon spectrum by summing both the hard perturbative QCD photons together with the thermal emission. We also compare to the $`e^+e^{}`$ spectrum as measured at the highest multiplicities. The background cocktail is the same as used by the experiment itself. The thermal rates are those computed by Eletsky, Belkacem, Ellis and Kapusta  based on kinetic theory and hadronic data. These rates are nearly identical to those computed by Rapp and Wambach ; both of them show a significant shift in spectral weight from the rho-meson region to lower masses. The results are shown in the following ten figures. Further details may be found in ref. .

## 2 Conclusions

The space-time evolution of coarse–grained UrQMD is very similar to that of hydrodynamics. Both represent the hadronic data quite well. The thermal photon and dilepton rates are relatively well understood. The theoretical predictions for both photons and dileptons agree rather well with data for central heavy ion collisions at 158 GeV per nucleon, although there may be a slight underestimate of dileptons in the mass range from 500 to 700 MeV. We should celebrate the success of both theory and experiment!


text:
# Mixing Time of the Rudvalis Shuffle

## 1. Introduction

In earlier work (Wilson, 2001) we derived upper and lower bounds on the mixing time of a variety of Markov chains, including Markov chains on lozenge tilings, card shuffling, and exclusion processes. The mixing time of a Markov chain is the time it takes to approach its stationary distribution, which is often measured in total variation distance (defined below). In this article we focus on the method for lower bounding the mixing time, and extend its applicability to the Rudvalis card shuffling Markov chain (defined below) and related shuffles.

Let $`P_x^t`$ denote the distribution of the Markov chain started in state $`x`$ after it is run for $`t`$ steps, and let $`\mu `$ denote the stationary distribution of the Markov chain. The total variation distance between distributions $`P^t`$ and $`\mu `$ is defined by

$$P_x^t\mu _{\mathrm{TV}}=\underset{A}{\mathrm{max}}\left|P_x^t(A)\mu (A)\right|=\frac{1}{2}\underset{y}{}\left|P_x^t(y)\mu (y)\right|=\frac{1}{2}P^t\mu _1,$$

and the mixing time is the time it takes for $`\mathrm{max}_xP_x^t\mu _{\mathrm{TV}}`$ to become small, say smaller than $`\epsilon `$. See (Aldous and Fill, 2005; Diaconis, 1988, 1996) for further background.

Arunas Rudvalis proposed the following shuffle: with probability $`1/2`$ move the top card to the bottom of the deck, and with probability $`1/2`$ move it to the second position from the bottom. Hildebrand (1990) showed that the Rudvalis shuffle mixes in $`O(n^3\mathrm{log}n)`$ time. Diaconis and Saloff-Coste (1995) studied a variation, the shift-or-swap shuffle, which at each step either moves the top card to the bottom of the deck or exchanges the top two cards, each move with probability $`1/2`$. Diaconis and Saloff-Coste (1993) also studied a symmetrized version of the Rudvalis shuffle, which at each step does one of four moves each with probability $`1/4`$: move top card to bottom, move bottom card to top, exchange top two cards, or do nothing. In each case a $`O(n^3\mathrm{log}n)`$ upper bound on the mixing time is known, but order $`n^3\mathrm{log}n`$ lower bounds were not known.

To lower bound the mixing time, one finds a set $`A`$ of states such that $`P^t(A)`$ is close to $`1`$ and $`\mu (A)`$ is close to $`0`$. The approach taken in (Wilson, 2001) uses an eigenvector $`\mathrm{\Phi }`$ of the Markov chain. If $`X_t`$ denotes the state of the Markov chain at time $`t`$, then $`𝔼[\mathrm{\Phi }(X_{t+1})X_t]=\lambda \mathrm{\Phi }(X_t)`$. To obtain a good lower bound, we need $`\lambda <1`$ but $`\lambda 1`$. Since $`\lambda <1`$, in stationarity $`𝔼[\mathrm{\Phi }(X)]=0`$, but since $`\lambda 1`$, it takes a long time before $`𝔼[\mathrm{\Phi }(X_t)]0`$. If furthermore the eigenvector is “smooth” in the sense that $`𝔼[|\mathrm{\Phi }(X_{t+1})\mathrm{\Phi }(X_t)|^2X_t]`$ is never large, then we can bound the variance $`\mathrm{\Phi }(X_t)`$, showing that it is with high probability confined to a small interval about its expected value. Then provided that $`𝔼[\mathrm{\Phi }(X_t)]`$ is large enough, we can reliably distinguish $`\mathrm{\Phi }(X_t)`$ from $`\mathrm{\Phi }(X)`$ in stationarity, which implies that the Markov chain has not yet mixed by time $`t`$. Saloff-Coste (2002) gives an exposition of this and related ideas.

For the Rudvalis card shuffling Markov chain and its variants, there are a few difficulties when directly applying this approach to lower bound the mixing time. The eigenvectors that one wants to use are complex-valued rather than real-valued, and $`\mathrm{\Phi }(X_t)`$ is no longer confined to a small interval around $`𝔼[\mathrm{\Phi }(X_t)]`$. Instead what happens is that $`\mathrm{\Phi }(X_t)`$ is with high probability confined to a narrow annulus centered at $`0`$. $`𝔼[\mathrm{\Phi }(X_t)]`$ becomes too small too quickly, and $`\mathrm{Var}[\mathrm{\Phi }(X_t)]`$ remains too large to be useful.

To lower bound the mixing time we want to in effect work with $`|\mathrm{\Phi }(X_t)|`$ and forget about $`\mathrm{arg}\mathrm{\Phi }(X_t)`$. To do this we start by lifting the Markov chain to a larger state space, and let us denote the state at time $`t`$ of the lifted chain by $`(X_t,Y_t)`$. (The mixing time of the lifted Markov chain will upper bound the mixing time of the original chain, so at the outset it is not clear that we can lower bound the mixing time of the original chain by considering its lifted version.) We find an eigenvector $`\mathrm{\Psi }`$ on the lifted chain such that for all $`x`$, $`y_1`$ and $`y_2`$, $`|\mathrm{\Psi }(x,y_1)|=|\mathrm{\Psi }(x,y_2)|`$, so that $`|\mathrm{\Psi }(X_t)|`$ is well-defined. If we show that $`\mathrm{\Psi }(X_t,Y_t)`$ is with high probability close to $`𝔼[\mathrm{\Psi }(X_t,Y_t)]`$, which in turn is far from $`0`$, it will follow that $`|\mathrm{\Psi }(X_t)|`$ is with high probability confined to a small interval far from $`0`$, making it statistically distinguishable from $`|\mathrm{\Psi }(X)|`$ in stationarity, implying that the Markov chain has not mixed by time $`t`$.

In the following sections we carry out these ideas to obtain lower bounds on the mixing time that match (to within constants) the previously obtained upper bounds. Specifically, we show

###### Theorem 1.

For any fixed $`\epsilon >0`$, after $`\frac{1o(1)}{8\pi ^2}n^3\mathrm{log}n`$ shuffles of the Rudvalis shuffle, $`\frac{1o(1)}{2\pi ^2}n^3\mathrm{log}n`$ shuffles of the shift-or-swap shuffle, or $`\frac{1o(1)}{\pi ^2}n^3\mathrm{log}n`$ shuffles for the symmetrized Rudvalis shuffle, the distribution of the state of the deck has variation distance $`1\epsilon `$ from uniformity.

## 2. Lifting the Shuffles

When the top card is placed at the bottom of the deck, the position of any given card is cycically shifted left, so we will call this move “shift-left”, and similarly “shift-right” is the move which places the bottom card on top of the deck. The move which exchanges the top and bottom cards will be called “swap”, and the move which does nothing will be called “hold”. Thus the moves of the Rudvalis Markov chain are “shift-left” and “swap & shift-left”, while the moves of the variation considered by Diaconis and Saloff-Coste are “shift-left” and “swap”, and the moves of the symmetrized version are “shift-left”, “shift-right”, “swap”, and “hold”.

The state $`X_t`$ of the Markov chain is the permutation giving the order of the cards at time $`t`$. Let $`\mathrm{}`$ (the diamond-suit symbol) denote a particular card of interest, and let $`X_t(\mathrm{})`$ denote the location of card $`\mathrm{}`$ within the deck at time $`t`$, where the positions are numbered from $`1`$ to $`n`$ starting from the top of the deck. When the Markov chain does a shift-left or shift-right, while the position of a card $`\mathrm{}`$ will change, all the cards get moved together, so it does not have such a large randomizing effect on the permutation. We will track the position of card $`\mathrm{}`$, but we should also track the amount of shifting. So when we lift the Markov chain to $`(X_t,Y_t)`$, the lifted Markov chain will also keeps track of

$$Y_t=\text{\# shift-left’s}\text{\# shift-right’s}modn.$$

For the Rudvalis shuffle $`Y_t=tmodn`$ deterministically, whereas for the other two variations, $`Y_t`$ will be a random number between $`1`$ and $`n`$ which approaches uniformity in $`O(n^2)`$ time.

Recall that we need an eigenvector $`\mathrm{\Psi }`$ of the lifted chain $`(X_t,Y_t)`$ such that $`|\mathrm{\Psi }(X_t,Y_t)|`$ is a function of $`X_t`$ alone. For a given card $`\mathrm{}`$, let

$$\mathrm{\Psi }_{\mathrm{}}(X_t,Y_t)=v(X_t(\mathrm{}))\mathrm{exp}(Z_t(\mathrm{})2\pi i/n),$$

where

$$Z_t(\mathrm{})=X_t(\mathrm{})X_0(\mathrm{})+Y_tmodn,$$

and $`v()`$ is a function, to be determined later, which makes $`\mathrm{\Psi }_{\mathrm{}}`$ an eigenvector. Initially $`Z_0(\mathrm{})=0`$, and the only time that the $`Z_t(\mathrm{})`$ changes is when the card $`\mathrm{}`$ gets transposed. The dynamics of $`(X_t(\mathrm{}),Z_t(\mathrm{}))`$ (mod $`n`$) are summarized by

$$(X_{t+1}(\mathrm{}),Z_{t+1}(\mathrm{}))=\{\begin{array}{cc}(X_t(\mathrm{}),Z_t(\mathrm{}))\hfill & \text{if move was “hold”}\hfill \\ (X_t(\mathrm{})1,Z_t(\mathrm{}))\hfill & \text{if move was “shift-left”}\hfill \\ (X_t(\mathrm{})+1,Z_t(\mathrm{}))\hfill & \text{if move was “shift-right”}\hfill \\ (X_t(\mathrm{})1,Z_t(\mathrm{})1)\hfill & \text{if move was “swap” and }X_t(\mathrm{})=1\hfill \\ (X_t(\mathrm{})+1,Z_t(\mathrm{})+1)\hfill & \text{if move was “swap” and }X_t(\mathrm{})=n\hfill \\ (X_t(\mathrm{}),Z_t(\mathrm{}))\hfill & \text{if move was “swap” and }\mathrm{}\text{ elsewhere}\hfill \end{array}$$

We define

$$\mathrm{\Psi }(X_t,Y_t)=\underset{\mathrm{}=1}{\overset{n}{}}\mathrm{\Psi }_{\mathrm{}}(X_t,Y_t).$$

If we increment $`y`$ while holding $`x`$ fixed, then $`\mathrm{\Psi }(x,y)`$ gets multiplied by the phase factor $`\mathrm{exp}(2\pi i/n)`$, so we have an eigenvector satisfying our requirement that $`|\mathrm{\Psi }(X_t,Y_t)|`$ be a function of $`X_t`$ alone.

## 3. The Lower Bound Lemma

The lower bounding lemma that we shall use is similar to Lemma 4 of (Wilson, 2001), but with the modifications described in the introduction. Saloff-Coste (2002) also gives a generalization of Lemma 4 from (Wilson, 2001) that may be used when the eigenvalues are complex, but the extension below seems to be better suited for the shuffles considered here.

###### Lemma 2.

Suppose that a Markov chain $`X_t`$ has a lifting $`(X_t,Y_t)`$, and that $`\mathrm{\Psi }`$ is an eigenfunction of the lifted Markov chain: $`𝔼[\mathrm{\Psi }(X_{t+1},Y_{t+1})(X_t,Y_t)]=\lambda \mathrm{\Psi }(X_t,Y_t)`$. Suppose that $`|\mathrm{\Psi }(x,y)|`$ is a function of $`x`$ alone, $`|\lambda |<1`$, $`\mathrm{}(\lambda )1/2`$, and that we have an upper bound $`R`$ on $`𝔼[|\mathrm{\Psi }(X_{t+1},Y_{t+1})\mathrm{\Psi }(X_t,Y_t)|^2(X_t,Y_t)]`$. Let $`\gamma =1\mathrm{}(\lambda )`$. Then when the number of Markov chain steps $`t`$ is bounded by

$$t\frac{\mathrm{log}\mathrm{\Psi }_{\mathrm{max}}+\frac{1}{2}\mathrm{log}\frac{\gamma \epsilon }{4R}}{\mathrm{log}(1\gamma )},$$

the variation distance of $`X_t`$ (the state of the original Markov chain) from stationarity is at least $`1\epsilon `$.

The proof of this modified lemma is similar to the proof of Lemma 4 in (Wilson, 2001), but for the reader’s convenience we give the modified proof. In the following sections we determine the functions $`v()`$ for the Markov chains which give the requisite eigenfunction $`\mathrm{\Psi }`$, and then use Lemma 2 to obtain the mixing time bounds stated in Theorem 1.

###### Proof of Lemma 2.

Let $`\mathrm{\Psi }_t=\mathrm{\Psi }(X_t,Y_t)`$, and $`\mathrm{\Delta }\mathrm{\Psi }=\mathrm{\Psi }_{t+1}\mathrm{\Psi }_t`$. By induction

$$𝔼[\mathrm{\Psi }_t(X_0,Y_0)]=\mathrm{\Psi }_0\lambda ^t.$$

By our assumptions on $`\lambda `$, in equilibrium $`𝔼[\mathrm{\Psi }]=0`$.

We have $`𝔼[\mathrm{\Delta }\mathrm{\Psi }(X_t,Y_t)]=(\lambda 1)\mathrm{\Psi }_t`$ and

$`\mathrm{\Psi }_{t+1}\mathrm{\Psi }_{t+1}^{}`$ $`=\mathrm{\Psi }_t\mathrm{\Psi }_t^{}+\mathrm{\Psi }_t\mathrm{\Delta }\mathrm{\Psi }^{}+\mathrm{\Psi }_t^{}\mathrm{\Delta }\mathrm{\Psi }+|\mathrm{\Delta }\mathrm{\Psi }|^2`$
$`𝔼[\mathrm{\Psi }_{t+1}\mathrm{\Psi }_{t+1}^{}(X_t,Y_t)]`$ $`=\mathrm{\Psi }_t\mathrm{\Psi }_t^{}[1+(\lambda 1)^{}+(\lambda 1)]+𝔼[|\mathrm{\Delta }\mathrm{\Psi }|^2X_t]`$
$`\mathrm{\Psi }_t\mathrm{\Psi }_t^{}[2\mathrm{}(\lambda )1]+R`$
and so by induction,
$`𝔼[\mathrm{\Psi }_t\mathrm{\Psi }_t^{}]`$ $`\mathrm{\Psi }_0\mathrm{\Psi }_0^{}[2\mathrm{}(\lambda )1]^t+{\displaystyle \frac{R}{22\mathrm{}(\lambda )}},`$
then subtracting $`𝔼[\mathrm{\Psi }_t]𝔼[\mathrm{\Psi }_t]^{}`$,
$`\mathrm{Var}[\mathrm{\Psi }_t]`$ $`\mathrm{\Psi }_0\mathrm{\Psi }_0^{}\left[[2\mathrm{}(\lambda )1]^t(\lambda \lambda ^{})^t\right]+{\displaystyle \frac{R}{22\mathrm{}(\lambda )}}.`$
Since $`(1\lambda )(1\lambda ^{})0`$, we have $`\lambda \lambda ^{}2\mathrm{}(\lambda )1`$, and by assumption $`2\mathrm{}(\lambda )10`$. Hence $`(\lambda \lambda ^{})^t[2\mathrm{}(\lambda )1]^t`$, and we have for each $`t`$
$`\mathrm{Var}[\mathrm{\Psi }_t]`$ $`{\displaystyle \frac{R}{22\mathrm{}(\lambda )}}={\displaystyle \frac{R}{2\gamma }}.`$
From Chebychev’s inequality,

$$\mathrm{Pr}\left[|\mathrm{\Psi }_t𝔼[\mathrm{\Psi }_t]|\sqrt{R/(2\gamma \epsilon )}\right]\epsilon .$$

As $`𝔼[\mathrm{\Psi }_{\mathrm{}}]=0`$, if $`𝔼[\mathrm{\Psi }_t]\sqrt{4R/(\gamma \epsilon )}`$, then the probability that $`|\mathrm{\Psi }_t|`$ deviates below $`\sqrt{R/(\gamma \epsilon )}`$ is at most $`\epsilon /2`$, and the probability that $`|\mathrm{\Psi }|`$ in stationarity deviates above this threshold is at most $`\epsilon /2`$, so the variation distance between the distribution at time $`t`$ and stationarity must be at least $`1\epsilon `$. If we take the initial state to be the one maximizing $`\mathrm{\Psi }_0`$, then

$$𝔼[|\mathrm{\Psi }_t|]=|\mathrm{\Psi }_{\mathrm{max}}||\lambda |^t|\mathrm{\Psi }_{\mathrm{max}}|(\mathrm{}(\lambda ))^t=|\mathrm{\Psi }_{\mathrm{max}}|(1\gamma )^t\sqrt{4R/(\gamma \epsilon )}$$

when

$`t`$ $`{\displaystyle \frac{\mathrm{log}\left[\mathrm{\Psi }_{\mathrm{max}}÷\sqrt{\frac{4R}{\gamma \epsilon }}\right]}{\mathrm{log}(1\gamma )}}.\mathit{}`$
## 4. The Rudvalis Shuffle

The first shuffle we consider is the original Rudvalis Markov chain. It will be instructive to consider a slight generalization, where the swap & shift-left move takes place with probability $`p`$, and the shift-left move takes place with probability $`1p`$. We shall assume that $`0<p<1`$ and that $`p`$ is independent of $`n`$. The particular values of $`p`$ that we are interested in are $`p=1/2`$ (for the original Rudvalis chain) and $`p=1/3`$.

We need to find an eigenvector for the random walk that a single card takes under this shuffle. We remark that this random walk is similar in nature to (but distinct from) a class of random walks, known as daisy chains, for which Wilmer (1999) obtained eigenvalues and eigenvectors. From other work of Wilmer (2002), it readily follows that the position of a single card takes order $`n^3`$ steps to randomize, and that the precise asymptotic distance from stationarity of the card’s position after $`cn^3`$ shuffles is given by an explicit expression involving theta functions.

###### Lemma 3.

The random walk followed by a card $`\mathrm{}`$ under the lifted Rudvalis shuffle has an eigenvector of the form

$$\mathrm{\Psi }_{\mathrm{}}(x,z)=v(x)e^{2\pi iz/n}$$

where $`v(x)`$ is the $`x^{\text{th}}`$ number in the list

$$\lambda ^{n2},\mathrm{},\lambda ^2,\lambda ,1,\chi ,$$

the eigenvalue is

$$\lambda =1\frac{p}{1p}\frac{4\pi ^2}{n^3}+O(1/n^4),$$

and

$$\chi =1+\frac{p}{1p}\frac{2\pi i}{n}+O(1/n^2).$$

###### Proof.

Let $`w=\mathrm{exp}(2\pi i/n)`$. If at time $`t`$ card $`\mathrm{}`$ is in any location between $`2`$ and $`n1`$, then

$$\mathrm{\Psi }_{\mathrm{}}(X_{t+1},Y_{t+1})=\lambda \mathrm{\Psi }_{\mathrm{}}(X_t,Y_t)$$

deterministically. To ensure that

$$𝔼[\mathrm{\Psi }_{\mathrm{}}(X_{t+1},Y_{t+1})(X_t,Y_t)]=\lambda \mathrm{\Psi }_{\mathrm{}}(X_t,Y_t)$$

when $`X_t(\mathrm{})=1`$, we require

$`pw^1+(1p)\chi `$ $`=\lambda ^{n1}`$
$`\chi `$ $`={\displaystyle \frac{\lambda ^{n1}pw^1}{1p}},`$
and for when $`X_t(\mathrm{})=n`$ we need

$`pw\chi +(1p)`$ $`=\lambda \chi `$
$`\chi `$ $`={\displaystyle \frac{1p}{\lambda pw}}.`$
Given these two equations, $`\mathrm{\Psi }_{\mathrm{}}`$ will be an eigenvector with eigenvalue $`\lambda `$. Thus,

$$f(\lambda )=\lambda ^npw\lambda ^{n1}pw^1\lambda 1+2p=0.$$

To identify a root of this polynomial, we use Newton’s method: $`z_{k+1}=z_kf(z_k)/f^{}(z_k)`$, starting with $`z_0=1`$. By Taylor’s theorem,

$$|f(z_{k+1})|\frac{1}{2}\underset{0u1}{\mathrm{max}}|f^{\prime \prime }(uz_k+(1u)z_{k+1})|\times \left|\frac{f(z_k)}{f^{}(z_k)}\right|^2.$$

If $`|z1|1/n^2`$, then $`f^{}(z)=(1p)n+O(1)`$ and $`f^{\prime \prime }(z)=(1p)n^2+O(n)`$. Consequently, if $`|z_k1|1/n^2`$ and $`|z_{k+1}1|1/n^2`$, then

$$|f(z_{k+1})|\frac{1+O(1/n)}{2}\frac{1}{1p}|f(z_k)|^2.$$

Since $`f(z_0)=p(2ww^1)=p4\pi ^2/n^2+O(1/n^4)`$, for large enough $`n`$ we have by induction that $`|f(z_k)|(1p)(p/(1p)4\pi ^2/n^2)^{2^k}`$, $`|z_{k+1}z_k|(p/(1p)4\pi ^2/n^2)^{2^k}/(n+O(1))`$, and $`|z_{k+1}z_0|O(1/n^3)`$. Thus, for large enough $`n`$, the sequence $`z_0,z_1,z_2,\mathrm{}`$ converges to a point $`\lambda `$, which by continuity, must be a zero of $`f`$. Since $`z_1=1p/(1p)4\pi ^2/n^3+O(1/n^4)`$ and $`|\lambda z_1|=O(1/n^5)`$, we conclude that the polynomial $`f`$ has a root at

$`\lambda `$ $`=1{\displaystyle \frac{p}{1p}}{\displaystyle \frac{4\pi ^2}{n^3}}+O(1/n^4).\mathit{}`$
It is easy to check that $`\mathrm{\Psi }_{\mathrm{max}}=n+O(1/n)`$. Next we evaluate $`R`$ for this eigenvector.

$$\frac{\mathrm{\Psi }_{\mathrm{}}(X_{t+1},Y_{t+1})\mathrm{\Psi }_{\mathrm{}}(X_t,Y_t)}{w^{Z_t(\mathrm{})}}=\{\begin{array}{cc}(\lambda 1)\lambda ^{X_t(\mathrm{})}=O(1/n^3)\hfill & \text{if }2X_t(\mathrm{})n1\hfill \\ \chi \lambda ^{n2}=O(1/n)\hfill & \text{if }X_t(\mathrm{})=1\text{, shift-left}\hfill \\ w^1\lambda ^{n2}=O(1/n)\hfill & \text{if }X_t(\mathrm{})=1\text{, swap \& shift-left}\hfill \\ 1\chi =O(1/n)\hfill & \text{if }X_t(\mathrm{})=n\text{, shift-left}\hfill \\ w\chi \chi =O(1/n)\hfill & \text{if }X_t(\mathrm{})=n\text{, swap \& shift-left}\hfill \end{array}$$

Adding up these contributions over the various cards $`\mathrm{}`$, we find

$`|\mathrm{\Psi }(X_{t+1},Y_{t+1})\mathrm{\Psi }(X_t,Y_t)|`$ $`O(1/n)`$
$`R=𝔼[|\mathrm{\Psi }(X_{t+1},Y_{t+1})\mathrm{\Psi }(X_t,Y_t)|^2(X_t,Y_t)]`$ $`O(1/n^2).`$
Plugging $`\lambda `$, $`\mathrm{\Psi }_{\mathrm{max}}`$, and $`R`$ into the Lemma 2 gives, for fixed values of $`\epsilon `$, a mixing time lower bound of

$$(1o(1))\frac{1p}{p}\frac{1}{8\pi ^2}n^3\mathrm{log}n.$$

## 5. The Shift-or-Swap Shuffle

At this point there are two ways we can approach the shift-or-swap shuffle. We can either take a direct approach in the same manner as in the previous section, or we can do a comparison with the Rudvalis shuffle with $`p=1/3`$.

If we take the direct approach, then we let $`v(x)`$ denote the $`x^{\text{th}}`$ element of the list

$$(2\lambda 1)^{n2},\mathrm{},(2\lambda 1)^2,2\lambda 1,1,\chi .$$

The constraints on $`\chi `$ are

$$\chi =\frac{2\lambda }{1+w^1}(2\lambda 1)^{n2}$$

and

$$\chi =\frac{1+w(2\lambda 1)^{n2}}{2\lambda }.$$

As in section 4, we solve for $`\lambda `$ and find that $`\lambda =1(1+o(1))\pi ^2/n^3`$, compute $`\mathrm{\Psi }_{\mathrm{max}}=\mathrm{\Theta }(n)`$ and $`R=O(1/n^2)`$, and obtain the mixing time lower bound of $`\frac{1o(1)}{2\pi ^2}n^3\mathrm{log}n`$ shuffles.

Alternatively, we can couple the shift-or-swap shuffle with the Rudvalis shuffle. Whenever the shift-or-swap shuffle makes a shift, the number of swap’s since the previous shift will be odd with probability $`1/3`$. If it is odd, then this is equivalent to a swap-&-shift-left move, and if it is even, then it is equivalent to a shift-left move. This explains why we were interested in the case $`p=1/3`$ in the previous section. After $`t`$ steps, with high probability $`(1+o(1))t/2`$ shift moves occured, which means that the state of the deck is what it would be after $`(1+o(1))t/2`$ Rudvalis shuffles (with $`p=1/3`$), possibly with an extra swap move. The lower bound for the shift-or-swap shuffle does not follow from the lower bound itself for the Rudvalis shuffle, but it does follow from what we showed about $`|\mathrm{\Psi }|`$ for the Rudvalis shuffle.

## 6. Symmetrized Version of the Rudvalis Shuffle

When analyzing the symmetrized version of the Rudvalis shuffle, it will be convenient to have symmetric coordinates, so we re-index the card locations to run from $`(n1)/2`$ up to $`(n1)/2`$, and the swaps occur at locations $`(n1)/2`$ and $`(n1)/2`$.

###### Lemma 4.

The random walk followed by a card $`\mathrm{}`$ under the lifted symmetrized Rudvalis shuffle has an eigenvector of the form

$$\mathrm{\Psi }_{\mathrm{}}(x,z)=v(x)e^{2\pi iz/n}$$

where

$$v(x)=\frac{1+\delta }{2}e^{i\theta x}+\frac{1\delta }{2}e^{i\theta x}=\mathrm{cos}(\theta x)+i\delta \mathrm{sin}(\theta x),$$

$$\theta =(1+o(1))\sqrt{2}\pi n^{3/2},$$

$$\delta =(1+o(1))\frac{1}{\sqrt{2}n^{1/2}},$$

and both $`\delta `$ and $`\theta `$ are real. The eigenvalue is

$$\lambda =\frac{1+\mathrm{cos}\theta }{2}=1\frac{\pi ^2+o(1)}{2}n^3.$$

###### Proof.

When $`x\pm (n1)/2`$, we can readily compute the eigenvalue $`\lambda `$ to be

$`\lambda `$ $`={\displaystyle \frac{\frac{1}{4}v(x+1)+\frac{1}{2}v(x)+\frac{1}{4}v(x1)}{v(x)}}`$
$`={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{cos}(\theta x)\mathrm{cos}\theta +i\delta \mathrm{sin}(\theta x)\mathrm{cos}\theta }{\mathrm{cos}(\theta x)+i\delta \mathrm{sin}(\theta x)}}+{\displaystyle \frac{1}{2}}`$
$`={\displaystyle \frac{1+\mathrm{cos}\theta }{2}}.`$
In order for our guessed eigenvector to be correct, there is also a constraint at $`x=(n1)/2`$:

$`{\displaystyle \frac{v(\frac{n1}{2})}{4}}+{\displaystyle \frac{v(\frac{n3}{2})}{4}}+(1+w){\displaystyle \frac{v(\frac{n1}{2})}{4}}`$ $`=\lambda ={\displaystyle \frac{v(\frac{n1}{2})}{2}}+{\displaystyle \frac{v(\frac{n3}{2})}{4}}+{\displaystyle \frac{v(\frac{n+1}{2})}{4}}`$
$`(1+w)v((n1)/2)`$ $`=v((n1)/2)+v((n+1)/2)`$
$`(1+w)(1+\delta )e^{i\theta (n1)/2}+(1+w)(1\delta )e^{i\theta (n1)/2}`$ $`=\begin{array}{cc}& (1+\delta )e^{i\theta (n1)/2}+(1\delta )e^{i\theta (n1)/2}+\hfill \\ & (1+\delta )e^{i\theta (n+1)/2}+(1\delta )e^{i\theta (n+1)/2}\hfill \end{array}`$
$`(w+2\delta +w\delta )e^{i\theta (n1)/2}+(w2\delta w\delta )e^{i\theta (n1)/2}`$ $`=(1+\delta )e^{i\theta (n+1)/2}+(1\delta )e^{i\theta (n+1)/2}`$
$`w\mathrm{cos}\frac{\theta (n1)}{2}\mathrm{cos}\frac{\theta (n+1)}{2}`$ $`=\delta \left[(2+w)i\mathrm{sin}\frac{\theta (n1)}{2}+i\mathrm{sin}\frac{\theta (n+1)}{2}\right].`$
The corresponding constraint at $`x=(n1)/2`$ is obtained by replacing $`w`$ with $`1/w`$ and replacing $`\theta `$ with $`\theta `$. Since these substitutions give the complex-conjugate of the above equation, the constraints at $`x=\pm (n1)/2`$ are equivalent.

Equating the real parts of this equation gives

$`\delta `$ $`={\displaystyle \frac{\mathrm{cos}\frac{2\pi }{n}\mathrm{cos}\frac{\theta (n1)}{2}\mathrm{cos}\frac{\theta (n+1)}{2}}{\mathrm{sin}\frac{2\pi }{n}\mathrm{sin}\frac{\theta (n1)}{2}}},`$
and equating the imaginary parts gives
$`\delta `$ $`={\displaystyle \frac{\mathrm{sin}\frac{2\pi }{n}\mathrm{cos}\frac{\theta (n1)}{2}}{2\mathrm{sin}\frac{\theta (n1)}{2}+\mathrm{cos}\frac{2\pi }{n}\mathrm{sin}\frac{\theta (n1)}{2}+\mathrm{sin}\frac{\theta (n+1)}{2}}}.`$
Cross-multiplying and performing trigonometric simplifications gives

$$\mathrm{sin}^2(\frac{2\pi }{n})\frac{\mathrm{sin}(\theta (n1))}{2}=\begin{array}{cc}& \mathrm{cos}^2(\frac{2\pi }{n})\frac{\mathrm{sin}(\theta (n1))}{2}\frac{\mathrm{sin}(\theta (n+1))}{2}+\mathrm{cos}(\frac{2\pi }{n})\mathrm{sin}\theta \hfill \\ & +\mathrm{cos}(\frac{2\pi }{n})\mathrm{sin}(\theta (n1))2\mathrm{sin}\frac{\theta (n1)}{2}\mathrm{cos}\frac{\theta (n+1)}{2}\hfill \end{array}$$

so that

(1)
$$\left(\frac{1}{2}+\mathrm{cos}\frac{2\pi }{n}\right)\mathrm{sin}(\theta (n1))\frac{1}{2}\mathrm{sin}(\theta (n+1))\mathrm{sin}(\theta n)+\left(1+\mathrm{cos}\frac{2\pi }{n}\right)\mathrm{sin}(\theta )=0.$$

Equation (1) is exact, but to estimate a solution, we perform a series expansion in $`\theta `$

$`0`$ $`=\begin{array}{cc}& \left[\left({\displaystyle \frac{1}{2}}+\mathrm{cos}{\displaystyle \frac{2\pi }{n}}\right)(n1)(n+1)/2n+1+\mathrm{cos}{\displaystyle \frac{2\pi }{n}}\right]\theta \hfill \\ & \left[\left({\displaystyle \frac{1}{2}}+\mathrm{cos}{\displaystyle \frac{2\pi }{n}}\right)(n1)^3(n+1)^3/2n^3+1+\mathrm{cos}{\displaystyle \frac{2\pi }{n}}\right]{\displaystyle \frac{\theta ^3}{6}}+O(n^4\theta ^5)\hfill \end{array}`$
$`0`$ $`=\left[{\displaystyle \frac{2\pi ^2}{n}}+O(1/n^2)\right]\theta +\left[6n^2+O(n)\right]{\displaystyle \frac{\theta ^3}{6}}+O(n^4\theta ^5).`$
While $`\theta =0`$ is a solution, our expression for $`\delta `$ has a singularity at $`\theta =0`$, so we seek a different solution. Ignoring the error terms suggests $`\theta \sqrt{2}\pi n^{3/2}`$. Since the function in (1) is real-valued, we can appeal to the intermediate value theorem to show that there is in fact a root at

$$\theta \sqrt{2}\pi n^{3/2}.$$

For this value of $`\theta `$ we have

$$\lambda =\frac{1+\mathrm{cos}\theta }{2}1\frac{\pi ^2}{2}n^3,$$

and (using the second equation for $`\delta `$)

$`\delta `$ $`{\displaystyle \frac{2\pi /n}{2n\theta /2+n\theta /2+n\theta /2}}{\displaystyle \frac{1}{\sqrt{2}n^{1/2}}}.\mathit{}`$
Again $`\mathrm{\Psi }_{\mathrm{max}}=(1+o(1))n`$. Next we estimate $`R`$. If there is a shift-left, then provided card $`\mathrm{}`$ is not in position $`(n1)/2`$, we have

$`\mathrm{\Delta }\mathrm{\Psi }_{\mathrm{}}`$ $`=(\mathrm{cos}\theta 1)[\mathrm{cos}(\theta x)+i\delta \mathrm{sin}(\theta x)]e^{2\pi iz/n}+\mathrm{sin}\theta [\mathrm{sin}(\theta x)i\delta \mathrm{cos}(\theta x)]e^{2\pi iz/n}`$
$`=O(\theta ^2)+O(\theta ^2x)+O(\theta \delta )=O(n^3)+O(n^2)+O(n^2)=O(n^2).`$
If card $`\mathrm{}`$ is in position $`(n1)/2`$, then

$$\mathrm{\Delta }\mathrm{\Psi }_{\mathrm{}}=2i\delta \mathrm{sin}(\theta (n1)/2)e^{2\pi iz/n}=O(n^1).$$

Adding up these contributions over the different cards, we find $`\mathrm{\Delta }\mathrm{\Psi }=O(n^1)`$. Likewise $`\mathrm{\Delta }\mathrm{\Psi }=O(n^1)`$ if the move was a shift-right. For transposes, $`\mathrm{\Delta }\mathrm{\Psi }_{\mathrm{}}`$ is nonzero for only two cards, and for these it is $`O(n^1)`$. Thus in all cases we have $`|\mathrm{\Delta }\mathrm{\Psi }|^2O(n^2)`$, and so $`RO(n^2)`$. Plugging our values of $`\lambda `$, $`\mathrm{\Psi }_{\mathrm{max}}`$, and $`R`$ into Lemma 2, we obtain, for fixed values of $`\epsilon `$, a lower bound on the mixing time of $`\frac{1o(1)}{\pi ^2}n^3\mathrm{log}n`$ shuffles.

## 7. Remarks

We have seen how to extend the lower bound technique used in (Wilson, 2001) to shuffles that are much slower than what the position of a single card would indicate. Interestingly, for the shift-or-swap and symmetrized-Rudvalis shuffles, the spectral gap for the lifted shuffle is smaller than the spectral gap of the shuffle itself, so it is curious that we obtained a lower bound for these shuffles by considering their lifted versions.

In an early draft, we lower bounded the mixing time of the original Rudvalis shuffle without considering its lifted version, and this earlier approach might be considered simpler. But it is not clear how to lower bound the symmetrized Rudvalis shuffle without lifting it, and our current approach has the advantage that the analyses for all three shuffles treated here are similar. The original Rudvalis shuffle and its lifting are isomorphic, and the earlier analysis is effectively a special case of the present analysis where the lifting is not explicit.

We suspect that the constants in the lower bounds of Theorem 1 are tight. Rudvalis asked if his shuffle was the slowest shuffle evenly supported on two generators; the lower bounds given here suggest that the shift-or-swap shuffle (for odd $`n`$) is slower by a factor of $`4`$.

## Acknowledgements

The author thanks Laurent Saloff-Coste and Persi Diaconis for calling attention to these shuffles, and Laurent Saloff-Coste for comments on an earlier draft.