text:
# The light curve of the semiregular variable L2 Puppis: I. A recent dimming event from dust

## 1 Introduction

L2 Puppis (HR 2748; HIP 34922) is a bright nearby red giant with a pulsation period of about 140 d. Its spectral type of M5eIII and luminosity of 1500 L indicate that it is evolving towards the tip of the Asymptotic Giant Branch (AGB). Evidence for mass loss at a rate of $`3\times 10^7\mathrm{M}_{}\mathrm{yr}^1`$ (Jura et al., 2002) supports this. L2 Pup is possibly the nearest star in this evolutionary phase, at a Hipparcos distance of $`61\pm 5`$ pc. Among known long-period AGB stars, only R Doradus has a similar distance. At 12 microns, L2 Pup is among the 15 brightest sources in the IRAS point source catalogue.

L2 Pup is unusual in several respects. Firstly, it shows a high degree of optical polarization, with a variable wavelength dependence that implies a long timescale for the growth and dissipation of dust grains (of the order of a decade; Magalhaes et al. 1986). Secondly, CO measurements by Kerschbaum & Olofsson (1999) indicate a very low expansion velocity (about 2.5 km s-1), which led them to label L2 Pup as an extreme case, with one of the smallest expansion velocities ever measured for an AGB star. The slow wind from L2 Pup led Winters et al. (2000, 2002) to suggest that this star could represent their B-model, in which mass loss is driven entirely by pulsations, without any significant input from radiation pressure on dust grains. This has been further discussed by Jura et al. (2002), who modelled the mass loss and suggested that the pulsations may be non-radial.

Thirdly, as we report here, this star has shown a remarkable change in mean visual magnitude over the past century, and is currently undergoing a dramatic dimming. We present visual and infrared photometry which characterizes this behaviour, and argue that the most likely cause is the formation of dust along the line of sight. We also present the first 10-micron spectrum of L2 Pup, obtained with the Japanese IRTS satellite, which shows strong silicate emission.

## 2 Light curve

### 2.1 Visual observations

We have analysed visual observations of L2 Pup from the following sources: the Royal Astronomical Society of New Zealand (RASNZ; 12100 measurements by 110 observers, including 1100 by A. Jones), the Variable Star Observers League in Japan (VSOLJ; 792 measurements by 6 observers) and the Association Francaise des Observateurs d’Etoiles Variables (AFOEV; 389 measurements by 3 observers). Only data from individual observers contributing 30 or more observations were used and we did not attempt to correct for offsets between observers. The top panel of Fig. 1 shows the combined data, binned to 10-day averages.

The variation in mean magnitude is evident. Before 1945 the data are too patchy and/or noisy to reach strong conclusions. However, it is clear that the star faded significantly after 1960, recovering around 1975 but remaining fainter by about 0.5 mag compared to the 1950s. Recently, starting around 1994, a dramatic fading has occurred. Even at maximum, the star now remains fainter than magnitude 6, and it is no longer a naked-eye variable.

We have used wavelet analysis to search for period and amplitude changes. Wavelets have been used previously to study long-period variables (e.g., Szatmáry et al. 1996; Bedding et al. 1998; Kiss et al. 1999). We used the weighted wavelet Z-transform (WWZ; Foster 1996), developed specifically for unevenly sampled data. We experimented with different values for the parameter $`c`$, which defines the tradeoff between time resolution and frequency resolution Foster (1996), and settled on $`c=0.005`$ as a good compromise. More details of the application of the WWZ transform to long-period variables are given by Bedding et al. (1998).

The lower three panels of Figure 1 show the wavelet plots for L2 Pup, based on the light curve in the top panel. The second panel shows the WWZ transform, with the grey scale indicating the significance of each frequency as a function of time (see Bedding et al. 1998). Only a small range of frequencies is shown – there is no evidence for significant power outside this range. The third and fourth panels show, for each time bin, the semi-amplitude (in magnitudes) and period (in days) corresponding to the peak of the WWZ in the second panel.

It is clear from Fig. 1 that the period of L2 Pup has remained within a fairly narrow range (135–145 d) for the past 75 years. Such period jitter is common among semiregular and Mira variables. There is a slight indication of a lengthening of period before the onset of the recent dimming, from 136 to 142 days, but the change (at most 4%) is well within the range of normal Mira variables. The amplitude, on the other hand, has changed significantly over the years, but the changes are apparently uncorrelated with the period jitter.

L2 Pup is classified as an SRb variable, a class of stars with poorly defined periodicity. In fact, the period of this star is remarkably stable. The characteristics of the wavelet plot show that L2 Pup should be classified as SRa, a class closely related to Miras, often differing only in having smaller amplitudes.

### 2.2 Early observations

L2 Pup was discovered to be variable by Gould in 1872 (Cannon & Pickering, 1907). Roberts (1893) estimated a period of 137.2 d from the dates of nine maxima observed by himself, Gould and Williams during 1872–92. Cannon & Pickering (1907), in the Second Catalogue of Variable Stars, quoted a period of 140.15 d, which they attributed to Roberts. Cannon & Pickering (1909) listed dates of maxima going back to 1872, and also adopted a period of 140.15 d. It therefore seems that the period of L2 Pup has remained within a narrow range since its discovery.

The visual magnitude at maximum during the late 1800s was about 3.6 (Williams, 1897; Cannon & Pickering, 1909), consistent with modern (1950s) out-of-decline values. Müller & Hartwig (1922) listed the magnitude range since discovery as 3.4–4.6 at maximum and 5.8–6.2 at minimum. This suggests that, between 1872 and 1918, L2 Pup was never as faint as it currently is.

### 2.3 Near-infrared photometry

We have taken $`JHKL`$ infrared photometry from Whitelock et al. (2000), together with more recent observations obtained with the same system (see Table 1). The infrared observations are much less frequent than the visual ones and do not sample the pulsation cycle very well. Nevertheless, the mean infrared magnitudes can be estimated, and we have shown these in Fig. 2 by fitting a low-order polynomial to each waveband.

The dimming at $`V`$ is followed at all infrared bands, although with reduced amplitude. The $`JK`$ colour has reddened during the dimming.

The mean $`K`$ magnitude before the dimming was $`2.34`$ (Table 4 of Whitelock et al. 2000). More recently, the magnitude dropped to $`2.0`$ or fainter111Bedding & Zijlstra (1998) gave $`K=2.65`$, but this was a transcription error and should have read $`2.15`$.. The Hipparcos distance allows comparison with the $`K`$-band period–luminosity (P–L) relation for nearby Miras (Whitelock & Feast, 2000):

$$M_K=3.47\mathrm{log}P+0.84.$$
(1)
For L2 Pup, this relation predicts $`M_K=6.64`$, which compares relatively poorly to the observed value of $`M_K=6.16`$. The recent fading has worsened the agreement. Thus, L2 Pup appears to be one of a small fraction of stars located below the P–L relation. We will return to this in Sec. 3.3.

Table 2 shows the mean colours and magnitudes before and after the dimming. The pre-dimming NIR magnitudes are taken from Table 4 of Whitelock et al. (2000). The current values are derived from the fits in Fig. 2, and are somewhat more uncertain because of the poor sampling of the pulsation cycle.

## 3 Discussion of the dimming

What causes the large variations in mean magnitude of L2 Pup, especially the current dramatic fading? We consider three possibilities: changes in temperature, luminosity and extinction.

### 3.1 Change in stellar temperature and/or luminosity

In late M-type giants, the optical flux is very sensitive to temperature. Even a slight cooling would lead to a sharp drop in $`V`$ via the formation of temperature-sensitive molecules (TiO, VO) that absorb strongly at optical wavelengths (Reid & Goldston, 2002).

We first need to estimate the effective temperature of L2 Pup. In the absence of a measured angular diameter, this requires reference to published temperature scales. The temperature scales for M-giants differ greatly for stars with and without extended atmospheres. Assuming no extended atmosphere, the pre-dimming $`VK`$ colour index of L2 Pup implies $`T_{\mathrm{eff}}3200\mathrm{K}`$ (Bessell et al., 1998; Houdashelt et al., 2000). Models with this temperature also give a $`JK`$ colour consistent with the observed value. Jura et al. (2002) adopted a temperature of $`3400\mathrm{K}`$, also taken from photospheric models, which is also consistent with the scale of van Belle et al. (1996).

However, stars with extended atmospheres have much lower colour temperatures (Feast, 1996; van Belle et al., 2002). The following relation was given by Feast (1996) for Miras, based on interferometric angular diameters:

$$\mathrm{log}T_{\mathrm{eff}}=0.474(JK)_0+4.059.$$
(2)
For L2 Pup, the pre-dimming $`JK`$ colour index yields the much lower value of $`T_{\mathrm{eff}}=2775K`$. Since this star is Mira-like and shows evidence for mass loss, the atmosphere is likely to be extended and the lower temperature is probably a better description of its optical/NIR spectral energy distribution.

During the fading, $`JK`$ increased to 1.55, which corresponds to a temperature decline of 660 K on the scale of Feast (1996). The increase in $`VK`$ to 8.5, would correspond to a temperature decline of about 300 K in the Bessell et al. (1998) models. The bolometric correction, BCK, increased from 3.05 to 3.21 on the low-temperature scale, and from 3.05 to 3.15 on the high temperature scale (Whitelock et al., 2000; Bessell et al., 1998). The luminosity decline in either case is about a factor of 1.8. Using $`LR^2T_{\mathrm{eff}}^4`$, we conclude that the radius would have decreased by 10% on the high-T scale and increased by 20% on the low-$`T`$ scale.

Such large changes in radius can be ruled out because of the stable period. The maximum period change allowed by the observations is about 6% (see Fig. 1). Since the period of a pulsating star scales inversely with the square root of stellar density, we can infer that the radius cannot have changed by more than 4% (since $`PR^{3/2}`$). In fact, we can take the relative stability of the pulsation period over the past 130 years to be evidence that the density of the star has not changed significantly. The period variations that are observed in L2 Pup are similar to those seen in many semiregular variables and are most likely due to the influence on the pulsation driving by random convective excitation (Christensen-Dalsgaard et al., 2001). We discuss this in more detail in Paper II (Bedding et al., 2002).

The period is stable on the same timescales over which the mean magnitude is variable, which argues that radius changes cannot be part of the explanation of the dimming. If we assume a constant radius, we could still produce the observed change in $`m_{\mathrm{bol}}`$ if the temperature decreased by about $`400\mathrm{K}`$. Such a temperature change cannot be ruled out but would require fine tuning. Furthermore, such a large change in stellar luminosity over a time scale of a decade only occurs in AGB stars in the immediate aftermath of a thermal pulse. Finding any star in such a short-lived phase is very unlikely. In any case, the thermal-pulse scenario can also be ruled out because, during this phase, the radius strongly contracts and the period shortens (e.g. Vassiliadis & Wood 1993), which is not observed in L2 Pup.

It is therefore unlikely that either the radius or the luminosity of the star has changed. This points towards variable extinction as the cause of the dimming.

### 3.2 Extinction from circumstellar dust

The plausible cause of the fading in L2 Pup is extinction by circumstellar dust, which would not affect the amplitude or period of pulsation of the underlying star. Extinction laws for interstellar dust are well studied and predict less extinction as one moves to longer wavelengths. This is certainly the case for L2 Pup, as can be seen from Fig. 2: the decline is much less dramatic in $`JHKL`$ than in $`V`$.

Direct evidence for extinction comes from the present $`JK`$ and $`KL`$, which are among the reddest of the M-type variables in the sample of Whitelock et al. (2000). All stars with such red colours have $`K`$-band amplitudes larger than 0.7 mag, while the amplitude for L2 Pup is only 0.29 mag at $`K`$. In fact, both $`\mathrm{\Delta }K`$ and $`\delta Hp`$ (0.71) are small compared to M-type variables. This large discrepancy between colour and amplitude suggests that the red colours are partly caused by extinction.

L2 Pup is also anomalously red for its period. The average relation between colour and period for Miras is

$$JK=0.39\pm 0.15+(0.71\pm 0.06)\mathrm{log}P,$$
(3)
(Whitelock et al., 2000). For L2 Pup, this relation predicts $`JK=1.13`$, compared with $`1.31`$ pre-dimming and $`1.55`$ at present. The observed values suggest that significant extinction existed even before the present dimming. Assuming a linear relation between the $`JK`$ excess and the visual extinction, and assuming that the post-1995 dimming of 2 mag at $`V`$ and 0.24 in $`JK`$ is due to extinction, the pre-dimming $`JK`$ excess of 0.18 would correspond to a visual extinction of about $`1.5`$ mag at $`V`$.

This estimate suggests that the long-term $`V`$-band variations may be due to variable extinction. The brightest epoch, in the 1950s, could in this case give the best indication of the unreddened magnitude. However, while L2 Pup is Mira-like in some properties, it is certainly not a true Mira and so results based on applying Eq. 3 should be treated with some caution.

### 3.3 Extinction curve and the Mira P–L relation

Typical extinction constants for interstellar dust are listed in Bessell et al. (1998). For the observed post-1995 $`JK`$ reddening, they give $`\mathrm{\Delta }A_K=0.15`$. This is less than half the decline at $`K`$ shown in Fig. 2, suggesting this extinction curve is too steep. To investigate this quantitatively, we show in Fig. 3 the locus $`(M,V)`$ of the mean infrared magnitude versus the $`V`$ magnitude, for $`M=J,H,K,L`$. Also shown are straight lines with slopes taken from Table 3 of Cardelli et al. (1989), which are standard values for extinction by interstellar dust. We can see that the observed slopes are in all cases larger than expected, and the discrepancy is progressively greater at longer wavelengths. The interstellar extinction curve appears inappropriate: instead the extinction appears to be ‘greyer.’ This could be due to a contribution from oxides, which can form at high temperatures (Henning, 1996). The infrared colour-colour diagram, shown in Fig. 4, confirms that the reddening is less steep than expected from interstellar extinction.

The colour differential suggests that $`A_\mathrm{K}/A_\mathrm{V}0.3`$. If we assume that the pre-dimming (1995) extinction at $`V`$ was $`A_\mathrm{V}1.5\mathrm{mag}`$, the $`K`$-band would already have been affected by 0.5 mag, so that $`M_{\mathrm{K},0}6.7`$. This is close to the value expected from the Mira P–L relation, and well within the 1-$`\sigma `$ uncertainty contributed by the error on the Hipparcos parallax. It therefore appears possible that the Mira P–L relation is an accurate indicator of the de-reddened magnitude of L2 Pup.

### 3.4 Dust location

For a pulsation velocity of 3 km s-1, the gas travels at most $`1R_{}`$ per year. The short timescale of the dimming thus suggests that the new dust is located within a few stellar radii of the star, possibly less. The rapid variation in polarization angle (Magalhaes et al., 1986) also suggests that dust forms close to the star. The timescale for dust formation derived by these authors based on the polarization (about a decade) is consistent with the timescale of extinction variations.

The polarization shows that the dust is not distributed isotropically, and the variation in polarization angle shows that the distribution is not constant. Dust may form continuously, but in different directions: the present dimming event would, in that case, correspond to dust formation along the line of sight.

### 3.5 Extinction episodes in other Mira variables

Dimming episodes lasting a number of pulsation cycles (decades) are seen in carbon Miras, with similarities to the dust obscuration events of R CrB stars (Percy et al., 1990). Lloyd Evans (1997) found that these episodes are accompanied by strong emission in the Swan C2 band 5165Å, while Whitelock et al. (1997) found them to occur in carbon stars with thick circumstellar envelopes.

Could the fading of the non-carbon star L2 Pup be a related phenomenon? Lloyd Evans (1997) listed five Miras with dust fadings which probably have relatively low C/O ratios (W Aql, U Cyg, Y Del, S Lyr and RU Vir), and Whitelock et al. (2000) added Y Vel to this list. Unlike L2 Pup, all these stars have long periods (433–490 d). And, in any case, inspection of the visual light curves of all six stars on the AAVSO Web site do not show any evidence for dimmings of the magnitude and longevity of those seen in L2 Pup.

Dimming episodes are seen in symbiotic Miras, where they are attributed to variable dust extinction (Whitelock, 1987). A typical dimming at $`J`$ is 1–2 mag, even larger than that seen in L2 Pup. We conclude that among non-symbiotic oxygen-rich Miras and semiregulars, the dimming of L2 Pup appears to be unique.

## 4 Mass loss

Whitelock et al. (2000) noted that L2 Pup stands out as one of only two low-amplitude SRs with a large $`K[12]`$ excess (the other being V CVn), implying a larger mass loss than is usual for SRs. There may be an effect from the low expansion velocity: the dust will remain close to the star longer, thus increasing the -$`\mu `$m excess.

Stronger evidence for on-going mass loss would come from a 10-$`\mu `$m silicate feature. Such a feature has not been reported, because of the lack of a published 10-micron spectrum. The lack of an original IRAS LRS spectrum is common for bright stars (e.g., $`o`$ Ceti). However, Volk & Cohen (1989) failed to recover the IRAS spectrum, which seems to have suffered from a technical glitch. ISO did not observe this star.

L2 Pup was included in the recent survey by the Japanese Infrared Telescope in Space (IRTS) satellite (Murakami et al., 1996). Observations were made using the Mid-Infrared Spectrometer (MIRS) (Roellig et al., 1994; Yamamura et al., 1996), which has an $`8\times 8`$ arcmin2 aperture and scanned with a speed of 4 arcmin/sec. Spectra were taken on 1995 March 30 and April 5, i.e., just before the dimming event. We re-performed the data reduction to correct for effects caused by the high brightness of the star. The coverage in this part of the sky by the IRTS was also below average, resulting in only two well-centred scans over L2 Pup. The averaged spectrum is shown in Fig. 5.

The spectrum shows a surprisingly strong silicate feature, with a band-to-continuum ratio of about 2.5. This give a silicate strength parameter (Schutte & Tielens, 1989) of $`A_{10}=2.5\mathrm{log}(F_{\mathrm{obs}}/F_{\mathrm{cont}})1`$. The continuum between 6 and 8 $`\mu `$m has a slope consistent with the stellar temperature but is brighter than expected from the $`L`$-band magnitude of $`2.9`$.

There is no evidence for a large dust halo around L2 Pup. We inspected a 4 square degree field in the four IRAS images made using the IRAS Software Telescope (Assendorp et al., 1995): the field shows complicated cirrus structure but no large envelope.

### 4.1 Mass-loss rate

The silicate feature is very strong compared to the limits on the underlying dust continuum. To quantify this, we fitted the spectrum using the dust model of Siebenmorgen & Kruegel (1992). The model uses as input a blackbody, which we assigned a temperature of $`T_{\mathrm{eff}}=2900\mathrm{K}`$ and luminosity $`L_{}=2.0\times 10^3\mathrm{L}_{}`$. The dust was described as a $`r^2`$ distribution (i.e., constant mass-loss rate) with an inner and outer radius: the mass-loss rate was converted to density assuming an expansion velocity $`v_{\mathrm{exp}}=2.5\mathrm{km}\mathrm{s}^1`$. The dust was included as silicate grains with a size distribution $`na^{3.5}`$, with diameter $`a>150`$ Å. In the fitting procedure, we adjusted the stellar temperature and luminosity, the dust density (or mass-loss rate) and inner radius, to reproduce the MIRS spectrum. The stellar parameters are largely determined by the shorter wavelength observations.

The 25 and 60-$`\mu `$m DIRBE measurements are much higher than the corresponding IRAS values. Part of this can be related to the different response curves: both these DIRBE pass bands cut off a few microns bluer than the IRAS bands. However, the increase over the predicted IRAS flux for the model spectrum is only about 10%. Although the DIRBE beam is much larger than that of IRAS (as the telescope was only 19cm, compared to 60cm for IRAS), but there are no other other bright point sources in the beam.

It is difficult to reconcile the photometry with the silicate feature. This is illustrated in Fig. 6, where the left panel shows a fit based the silicate feature, and the right panel shows a fit that includes all photometry. The two models differ in their outer radius, which is smaller for the right panel. The fitted parameters of the right panel (our preferred model) are listed in Table 3. The large number of free parameters means that the fits are not unique: i.e., a larger mass-loss rate with a thinner shell could also fit the spectrum.

The most constraining aspects of the spectrum are the lack of dust emission at 8$`\mu `$m and the strong silicate feature, which preclude the presence of hot dust. In our ‘best’ model, the inner radius is at a dust temperature of 350 K, well below the dust condensation temperature. The 17$`\mu `$m flux (also dominated by silicate emission) is also surprisingly sensitive to the inner radius and agrees with a dust shell detached from its condensation point. The outer radius is not as well determined as the inner radius. However, a single thin shell can fit the long-wavelength IRAS/DIRBE flux densities. There is no evidence for a cold distant dust shell, consistent with the low CO mass-loss rate found by Olofsson et al. (2002).

Between 4 and 8$`\mu `$m, the spectral energy distribution indicates a colour temperature consistent with the photospheric temperature. Jura et al. (2002) argued that dust emission dominates from 5 $`\mu `$m, based on evidence that the variability behaviour is very different longward of 5 $`\mu `$m. However, the emission in this band in Mira variables is dominated by water vapour at $``$2000 K, located in the extended envelope at about 1.5 stellar radii (Matsuura et al., 2002). This water vapour also raises the flux in this region above that extrapolated from the photospheric $`JHK`$ magnitudes, as seen in the fit. Our findings do not indicate a substantial contribution from dust to the spectrum shortward of 8 $`\mu `$m.

We derive a (pre-dimming) mass-loss rate of $`5\times 10^7\mathrm{M}_{}\mathrm{yr}^1`$. This is similar to that estimated by Jura et al. (2002) (post-dimming). They found a dust mass-loss rate 4 times larger than derived by us because they adopted a dust outflow velocity of 10 km s-1 rather than 2.5 km s-1, but used a much lower gas-to-dust ratio, which almost cancels this difference. Our high gas-to-dust ratio is based on the assumption that L2 Pup is a thick-disk star (see Sec. 4.3), for which we have assumed \[Fe/H\]$`=0.7`$. In reality the gas mass-loss rate is not well determined.

### 4.2 Silicate and Corundum

The model shows that the photometric data imply a narrower silicate feature than is observed. The most likely explanation is that the observed 10-$`\mu `$m feature is not purely silicate. In the scheme of Speck et al. (2000), this feature in L2 Pup would be classified either as ‘broad’ or possibly as ‘silicate D’. (The difference is subtle and requires continuum subtraction, for which we lack sufficient wavelength coverage.) They fitted the broad feature with a mixture of olivine (MgFeSiO4) and porous amorphous alumina (or corundum) (Al2O3). Corundum has the highest condensation temperature of any circumstellar solid, and is the first condensate to form at $`T1600`$ K. Silicates may form subsequently from gas-solid reactions with alumina (Tielens, 1990). Hron et al. (1997) suggested that the broad feature is seen in less evolved Mira variables. This may reflect the possiblity that, at low $`\dot{M}`$, the conversion of alumina to silicate (replacing the Al with Mg) remains incomplete, due to a low reaction rate.

We note that our model calculations used the pre-dimming spectrum. It would interesting to see whether the spectrum has changed significantly since the dimming. If alumina formed first, the peak of the 10$`\mu `$m feature may have shifted to the red. It is tempting to associate the alumina with the grey extinction curve, but it has little extinction in the optical (glass-like). Better candidates would be TiO2 or pyroxenes (Woitke, 1999). The precise agent responsible for the variable polarization also remains undetermined.

### 4.3 Evolutionary status

Jura et al. (2002) pointed out that the luminosity and effective temperature of L2 Pup do not unambiguously identify it as an AGB star, and allow the possibility that it is on the first ascent red giant branch. It is not clear whether the spectrum of L2 Pup shows technetium, which is usually taken as indicating an AGB star. Little et al. (1987) reported a probable detection, while Lebzelter & Hron (1999) listed only a possible detection.

However, if the $`K`$-band has significant obscuration then the bolometric magnitude has been underestimated. The bolometric correction would also have been overestimated because of the observed $`VK`$ is too large. The two effects combined would mean that the actual luminosity is substantially greater than has been assumed, and is perhaps sufficiently large that placement on the AGB is unambiguous.

Finally, as pointed out by Jura et al. (2002), the space motion of L2 Pup indicates it may be a thick disk star. Given that the thick disk accounts for up to 13% of the local star density (Chen et al., 2001), it is not unexpected to find a local Mira-like example. If L2 Pup is indeed a thick-disk star, as seems plausible but which cannot be stated with certainty, the progenitor would likely have had a low metallicity and a low mass.

## 5 Conclusions

Visual photometry of L2 Pup shows an unprecedented dimming over the past 5 years. The long-term light curve shows stable periodicity, and we argue that L2 Pup is Mira-like and should be classified as SRa. The period stability implies a constant stellar radius, which rules out temperature and/or luminosity variations as the cause of the dimming. Rather, the dimming seems to arise from an episode of dust formation close to the extended atmosphere. Episodic dust obscuration events are fairly common in carbon stars but have not been seen in (non-symbiotic) oxygen-rich stars. We suggest that dust forms continuously but anisotropically, with the current dimming event being due to to dust formation along the line of sight.

The red colours indicate reddening from dust, but the extinction curve is greyer than found for ISM dust. This could reflect a higher fraction of oxides. The change of colour during the dimming indicates that already before the dimming, the $`V`$-band magnitude was significantly affected by circumstellar or atmospheric extinction. L2 Pup was one of the few stars located below the Mira $`P`$$`L`$ relation: the derived $`K`$-band extinction puts the star in closer agreement with this relation.

We present a 10-$`\mu `$m spectrum showing strong silicate emission. These observations were carried out in 1995, just prior to the recent dimming. The silicate feature can be fitted with a detached, thin dust shell, with inner radius $`7\times 10^{14}\mathrm{cm}`$ and outer radius $`1.5\times 10^{15}\mathrm{cm}`$. Longer wavelength photometry shows no evidence for more distant, colder dust. We derive a mass-loss rate of $`\dot{M}_\mathrm{g}5\times 10^7\mathrm{M}_{}\mathrm{yr}^1`$, but this value depends on the assumed expansion velocity and metallicity — if the dust velocity is high, the actual mass-loss rate could be higher.

## Acknowledgments

We are extremely grateful to the observers who have provided visual observations, and to those who have maintained the databases and made them publicly available. We also thank Peter Tuthill for useful discussions. TRB and AR are grateful to the Australian Research Council for financial support.


format: text

text:
# A microscopic 2𝐷 lattice model of dimer granular compaction with friction

## I Introduction

The process of compaction in granular media attracts a great deal of attention in the scientific community, in particular amongst physicists and chemical engineers. In fact it displays features that are general enough to make it suitable for investigation through models that are relatively simple, at least for a non equilibrium process. Granular compaction therefore appears appealing for testing new promising and unifying ideas in the field of disordered systems .

Grains poured into a vessel fill it with loose arrangements and relatively low densities. The action of external perturbations, like shaking and tapping, in the presence of an external driving, like the gravity field, leads to a very slow increase of density through a rich phenomenology, displaying different regimes and both reversible and irreversible dynamical phases . A number of different lattice models have been recently investigated in order to unravel the microscopic mechanisms producing such behaviors, and to clarify to which extent they may be considered general, also in relation to other disordered systems like glasses .

In this work we introduce a lattice gas model for dimer compaction which includes, besides the vertical, also horizontal shaking and friction. Friction, although playing a fundamental role in real granular systems, has not been introduced so far in models for compaction. Horizontal shaking has been considered only experimentally for non Brownian spheres where it seems to favor local dense crystalline order . Highly anisotropic granulars, namely rods, have been recently investigated under vertical shaking. In contrast to sphere packings which tend to end up in disordered configurations, evidence is found that the particle anisotropy drives ordering . We try to incorporate in our model, albeit in a simplified manner, all these aspects, namely particle anisotropy, horizontal shaking and frictional dynamics, which we find to affect significantly the compaction process.

In particular, our model is based on the diffusional dynamics of dimers, considered as rigid and non overlapping particles which occupy two lattice sites on a square lattice. The model has no quenched disorder, but for this type of dynamics glass-like properties may be expected , because of the onset of geometrically frustrated configurations. Since the number of states with the closest packed density is exponentially large in the lattice size , this model is intrinsically different from the other lattice models recently investigated  and is suitable to understand how horizontal shaking and friction may influence the general characteristics attributed to the compaction process.

In the next section we describe in detail the model and the Monte Carlo method that we used. In Sec. III we discuss the results for the compaction dynamics without friction and illustrate the evolution of the packing structure. In Sec. IV we introduce the effect of friction and show its crucial role in the dynamical behavior. The final section is devoted to the conclusions.

## II Model

We consider a square lattice of $`N=L\times L`$ sites, with lateral periodic boundary conditions, an open boundary at the top and a rigid wall at the bottom. Elongated particles occupying two consecutive lattice sites (dimers) are inserted from the top one at the time with random (horizontal or vertical) orientation, letting them fall down keeping their orientation fixed, until reaching a stable position.

In this way we are able to prepare the system in an initial state which is saturated, i.e. no more particles can be put in, and with statistically reproducible density $`\rho _0`$, i.e. characterized by a precise mean value, $`\rho _00.587`$, corresponding to a random loose packing. In the present model, particles are subjected only to geometrical interactions, and the non-overlapping condition produces strong constraints on their relative positions (see Fig. 1). The aim of this paper is to study the system in presence of gravity and external vibrations, with the possibility of taking into account also a friction-like force for horizontal displacements. For this purpose we perform Monte Carlo (MC) simulations introducing a random diffusive dynamics in our model, which mimics the aforementioned forces and preserves the geometrical constraints. At this stage, we consider separately horizontal and vertical dimers and keep their orientation fixed without allowing rotations. We plan to introduce this feature, which would make the model more realistic, in the future. We consider in detail a single tap applied to the system for a fixed time $`\tau `$, comparing it with the multiple tapping for selected cases (see ). We first describe the dynamics without friction. The dynamics can be divided into two stages: $`(i)`$ for $`t<\tau `$ particles can move horizontally (left or right with probability $`p_h/2`$) and vertically (upwards with probability $`p_{up}`$ or downwards with probability $`p_{down}`$); $`(ii)`$ for $`t>\tau `$ only downward movement of particles is possible (i.e. $`p_{up}=p_h=0`$ and $`p_{down}=1`$). We note that $`\tau ^1`$ can be thought of as the quench rate of an initially annealed system. The probabilities for the different moves are normalized: $`p_h+p_{up}+p_{down}=1`$. Physically $`p_{down}`$ corresponds to the action of gravity, while $`p_h`$ and $`p_{up}`$ represent respectively the horizontal and vertical component of the external tap. The functional form chosen for the time dependence of $`p_h`$ and $`p_{up}`$ is

$$p_h(t)=p_h^0(1t/\tau )\theta (\tau t),$$
(1)
where $`\theta (x)`$ is the Heaviside function (the same expression holds for $`p_{up}(t)`$, with $`p_{up}^0`$ instead of $`p_h^0`$). $`p_{down}(t)`$ is determined from the normalization condition:

$$p_{down}(t)=1p_{up}(t)p_h(t)=p_{down}^0+(p_{up}^0+p_h^0)t/\tau $$
(2)
where $`p_{down}^01(p_{up}^0+p_h^0)`$. Since $`p_{down}(t)0`$ $`t`$ it must be $`p_{up}^0+p_{down}^01`$. If not stated explicitly otherwise, we assume $`p_{down}^0=p_{up}^0`$, so $`p_h^0`$ results to be $`p_h^0=12p_{up}^0`$, that is only one independent input parameter is needed. In each MC move one extracts a particle with uniform probability, chooses a move for this particle according to the values of $`p_h`$, $`p_{up}`$ and $`p_{down}`$, and performs the move if all the geometrical constraints are satisfied. One MC time step (MCS) corresponds to $`N`$ attempted moves. In the following, time is always given in units of MCS. To save CPU time we used an algorithm in which the attempted moves are always successful, and consistently updated time through probabilistic arguments (for the details of our computation see Ref.  on a reaction-diffusion model for dimers). We performed our MC simulations on a lattice with $`L=100`$, for which we checked that finite size effects are negligible.

As a next step we introduce friction between adjacent horizontal layers in the material in the following simplified way: the effective probability of making a horizontal move is set to

$$p_h^{eff}(t)=p_h(t)\mu (p_{down}(t)p_{up}(t))=$$

$$=[p_h^0\mu (1p_h^02p_{up}^0)[p_h^0+\mu (p_h^0+2p_{up}^0)]t/\tau ]\theta (t_0(\mu )t)$$
(3)
where

$$t_0(\mu )=\tau \left[\frac{p_h^0\mu (1p_h^02p_{up}^0)}{p_h^0+\mu (p_h^0+2p_{up}^0)}\right],$$
(4)
$`\mu `$ acts as a friction coefficient ($`\mu >0`$) and the “friction force” is proportional to the load, represented by the net vertical force, reducing the probability of the move. This assumption implies that the frictional force occurs only for horizontal moves, independently of whether an underlying dimer (i.e. a dimer in the adjacent lower row) is present. We are aware that this is a simplified approximation, but it is meant as a first step in describing this complex process. To ensure $`p_h^{eff}0`$ at $`t=0`$, $`\mu `$ has to satisfy the condition

$$\mu \frac{p_h^0}{(1p_h^02p_{up}^0)}.$$

In addition, for the friction force to give a negative contribution to $`p_h`$, one must have $`p_{down}(t)p_{up}(t)`$ $`t`$ (this is automatically satisfied with our choice $`p_{down}^0=p_{up}^0`$). Note that the normalization of the frictionless probabilities is still $`p_{up}+p_{down}+p_h=1`$. In other terms we have introduced a “sticking” probability (probability that the particle does not move) $`p_s`$, given by $`p_s=\mu (p_{down}p_{up})`$. If we put $`\mu =0`$ we recover the frictionless case.

## III Compaction dynamics without friction

We have analyzed the time behavior of the mean density of the system $`\rho (t)`$ (number of occupied sites normalized to the total number of sites) by measuring it in the lower $`30\%`$ of the box at time $`t`$. In this way we are sure to measure density at the bulk, since fluctuations can increase significantly in the proximity of the open boundary, resulting in a quite complicated compaction behavior in this region . Since the statistical fluctuations of the density can be quite relevant we have performed several MC realizations (up to $`500`$) of the process in order to obtain reliable results. After preparing the system in its low density configuration ($`\rho _0\rho (t=0)0.587`$), the diffusive dynamics starts as described in the previous section. $`p_h^0`$ is chosen as input parameter and consequently $`p_{up}^0=p_{down}^0=(1p_h^0)/2`$. We stop the simulations when a steady asymptotic value $`\rho _{\mathrm{}}`$ of the density is reached. We have studied the density evolution for different values of $`p_h^0`$ and of the shaking time $`\tau `$, in order to find out if some kind of scaling law, describing the behavior of the final density as a function of these parameters, exists. As far as we know, no systematic study concerning this point has been performed.

The time behavior of the density for several choices of $`\tau `$ is illustrated in Fig. 2, where a relatively low value of $`p_h^0`$ has been used ($`p_h^0=0.1`$).

We have observed that the dynamics gets drastically slower when $`\tau `$ increases; actually, because of our shaking rule, Eq. (1), when a long tap is applied the density does not significantly change in the first steps of the process, but the evolution takes place on a much wider time scale and finally a larger value of the asymptotic density is achieved. We have tried to fit different functional forms to the data in Fig. 2, looking in particular at the relaxation functions proposed in the seminal experimental paper of Ref. . It is clear that the observed dynamical behavior is very complex and is not compatible with a single relaxation time, i.e. a simple exponential. We have found that the commonly claimed inverse logarithmic relaxation does not hold for our system in the single tap case. Instead, the most suitable functional form for our data is a stretched exponential, as proposed by Nicodemi et al.  for the single tapping dynamics:

$$\rho (t;\tau ,p_h^0)=\rho _{\mathrm{}}(\tau ,p_h^0)C(\rho _{\mathrm{}}(\tau ,p_h^0)\rho _0)\mathrm{exp}[((t+t_0)/\tau _0)^\beta ]$$
(5)
where the fit parameters $`C`$, $`t_0`$, $`\tau _0`$ and $`\beta `$ depend in principle both on $`\tau `$ and $`p_h^0`$. Actually, it turns from our fits that $`C1.1`$ and $`\beta 10`$ are almost independent of $`\tau `$ and $`p_h^0`$. In Fig. 3(a) we show the density evolution for a small value of $`\tau `$.

As it can be seen, the intermediate-long time behavior can be accurately fitted by Eq. (5) (all fits have been performed for $`\rho 0.6`$). As $`\tau `$ increases, however, Eq. (5) fails to reproduce the long time regime. In fact, a power law tail develops, indicating that the compaction process slows down at long times for large shaking duration. The fitting function that we have used for the long time behavior is

$$\rho (t;\tau ,p_h^0)=\rho _{\mathrm{}}(\tau ,p_h^0)\frac{\rho _{\mathrm{}}(\tau ,p_h^0)\rho _0}{1+B(t/\tau _0^{})^\alpha }$$
(6)
This is illustrated in Fig. 3(b) for $`\tau =4000`$. Fit (6) has been applied for $`\tau >1000`$. The exponent $`\alpha `$ of this power law can be considered independent of $`\tau `$ ($`\alpha 5`$). However the stretched exponential function still describes the density behavior in the intermediate time regime well.

We have also studied the density relaxation for different values of $`p_h^0`$ at fixed $`\tau `$ ($`\tau =103`$). The corresponding plots are shown in Fig. 4, where it is clear that the role of the horizontal shaking amplitude is crucial in determining the asymptotic density.

Moreover we have considered the effect of different ratios $`p_{up}^0/p_{down}^0`$ on the density evolution, which we show in Fig. 5. When $`p_{up}^0>p_{down}^0`$ we observe a decompaction at the beginning of the process, which becomes more pronounced if the vertical tap is stronger.

This reflects the fact that, in the initial transient, density decreases at the bottom of the container for a strong tap. Such a decompaction upon vertical acceleration increase seems to be a genuine feature of two-dimensional systems, as claimed in . The saturation density is also sensitive to variations of $`p_{up}^0`$ and decreases for increasing $`p_{up}^0/p_{down}^0`$. However, the scaling behavior does not change and the density evolution can still be described by Eq. 5, but with a smaller value of $`\beta `$ with respect to the case $`p_{up}^0=p_{down}^0`$ ($`\beta 9.8`$ for $`p_{up}^0/p_{down}^0=1`$, $`\beta 7.6`$ for $`p_{up}^0/p_{down}^0=2`$ and $`\beta 6.5`$ for $`p_{up}^0/p_{down}^0=3`$).

In order to gain a better understanding of the degree of compaction of the system, we investigated the dependence of the saturation density $`\rho _{\mathrm{}}`$ on $`\tau `$ and $`p_h^0`$. The data for $`\rho _{\mathrm{}}`$ vs. $`\tau `$ are plotted in Fig. 6(a).

A crossover from a power law scaling to a logarithmic behavior is observed at a characteristic value $`\tau ^{}200`$:

$$\rho _{\mathrm{}}(\tau )=(\tau /\tau _l)^\delta \mathrm{for}\tau <\tau ^{}$$

$$\rho _{\mathrm{}}(\tau )=\rho _mB/\mathrm{ln}(\tau /\tau _r)\mathrm{for}\tau >\tau ^{}$$
(7)
where $`\tau _l4.510^3`$, $`\delta 0.073`$, $`\rho _m1`$, $`B0.6`$ and $`\tau _r18`$. This means that for large values of $`\tau `$ the process gets slower and slower, as observed before, eventually reaching a final value of the density $`\rho _m1`$ for $`\tau \mathrm{}`$ (which corresponds to about the maximal density for the system ). The variation of $`\rho _{\mathrm{}}`$ with respect to the ratio $`p_h^0/p_{down}^0`$ is illustrated in Fig. 6(b). The increase of $`\rho _{\mathrm{}}`$ is more pronounced for low values of $`p_h^0`$, while it tends to saturate afterwards. This compares qualitatively well to the experimental results for a horizontally shaken box filled with beads , where it was found that for low filling rates the packing crystallizes upon increase of the adimensional parameter $`\mathrm{\Gamma }=A\omega ^2/g`$ (where $`A`$ and $`\omega `$ are the amplitude and pulsation of the vibration and $`g`$ is the gravitational acceleration), which roughly corresponds to $`p_h^0/p_{down}^0`$.

In Fig. 7 we show some snapshots of one particular realization of the compaction process for $`\tau =4000`$.

We see how the packing evolves from a highly disordered state at $`t=0`$ to a polycrystalline state, made of many crystalline ordered domains, at the end of the compaction. As time increases the order proceeds from the bottom to the top of the packing while empty spaces reduce in size. The final state is therefore characterized by the presence of large clusters of horizontal and vertical dimers. In order to demonstrate that the shaking amplitude $`p_h^0`$ and time $`\tau `$ play a significant role, we have reported in Fig. 8 a comparison between the final configurations for $`p_h^0=0.1`$ and $`p_h^0=0.8`$ at fixed $`\tau =103`$ (upper part), and for $`\tau =80`$ and $`\tau =4000`$ at fixed $`p_h^0=0.1`$ (lower part).

We notice that an intense or/and long tap is effective in locally removing frozen configurations through collective rearrangements of particles, thus letting the defects migrate towards the free surface of the packing.

We have also investigated in detail the behavior of the parameters of the stretched exponential fit (Eq. (5)) $`t_0`$ and $`\tau _0`$, which is reported in Fig. 9.

We have observed a change of behavior in correspondence of the critical value of the shaking time $`\tau ^{}`$, the same value determining the change in the behavior of $`\rho _{\mathrm{}}`$ vs. $`\tau `$ described in Eq. (7). A similar transition was also signaled by Nicodemi et al. . The parameters of the stretched exponential relaxation Eq. (5) change according to:

$$t_0(\tau )=\tau /t_1+t_<\mathrm{for}\tau <\tau ^{}$$

$$t_0(\tau )=(\tau /t_2)^\gamma +a\mathrm{for}\tau >\tau ^{}$$
(8)
$$\tau _0(\tau )=\tau /\tau _1+\tau _<\mathrm{for}\tau <\tau ^{}$$

$$\tau _0(\tau )=(\tau /\tau _2)^\gamma +b\mathrm{for}\tau >\tau ^{}$$
(9)
where $`t_11.40`$, $`t_<130`$, $`t_2710^5`$, $`\gamma 0.47`$, $`a700`$, $`\tau _10.76`$, $`\tau _<150`$, $`\tau _210^4`$, $`b700`$. It is remarkable that both $`t_0`$ and $`\tau _0`$ vary linearly with $`\tau `$ for $`\tau <\tau ^{}`$ and both follow a power law with the same exponent $`\gamma `$ for $`\tau >\tau ^{}`$. This could mean that they are both determined by the same relaxation mechanism as $`\tau `$ changes (i.e. if we define $`\overline{t}_0=t_0t_<`$, $`\overline{\tau }_0=\tau _0\tau _<`$, $`\stackrel{~}{t}_0=t_0a`$ and $`\stackrel{~}{\tau }_0=\tau _0b`$, we see from Eqs. (8) and (9) that $`\overline{t}_0/\overline{\tau }_0`$ and $`\stackrel{~}{t}_0/\stackrel{~}{\tau }_0`$ are constants independent of $`\tau `$).

We note that the critical shaking time $`\tau ^{}`$ physically signals the crossover between two different regimes. In fact, a careful inspection of Fig. 2 shows that, for $`\tau <\tau ^{}`$, compaction occurs after the shaking (when $`p_h=p_{up}=0`$ and $`p_{down}=1`$) and is mainly driven by gravity. On the other hand, when $`\tau >\tau ^{}`$, the compaction time ($`t_{comp}`$), i.e. the time needed to reach the saturation density, is shorter than $`\tau `$, and vibrations act on the system till the end of the process (see Fig. 10).

We believe this is responsible for the power law tail for large values of $`\tau `$, since the combined action of shaking and gravity produces a slowing down of the dynamics, enabling at the same time to obtain denser packings. However, we are not sure whether this picture corresponds to a real dynamical phase transition in the system, and moreover the identification of an order parameter would be quite problematic. Therefore a deeper investigation on this point is needed and we hope to stimulate some future experimental and theoretical works.

Now we briefly mention the main changes occurring when a sequential tapping is applied to the system. Although the single tap dynamics can give insight on the dynamical evolution and the relaxation of the system, the sequential tapping procedure can be closely related to experimental situations. In Fig. 11(a) we compare the density evolution obtained by a single tap of duration $`\tau =10000`$ with that obtained by a sequential tapping applied for the same amount of time $`\tau `$. In particular we have generated $`5000`$ taps, each of which was applied for a time $`\tau /5000`$.

From the comparison we deduce that the multiple tapping dynamics is slower, suggesting that the fitting form Eq. 5 cannot describe properly the density evolution in this case. Besides, the saturation density is also lower. In Fig. 11(b) we characterize more quantitatively the dynamical behavior by plotting the density as a function of the number of taps. It turns out that in this case the intermediate-long time evolution of the density can be described by the claimed inverse logarithmic law, which was also found in the experiments :

$$\rho (t)=\rho _{\mathrm{}}c\frac{\rho _{\mathrm{}}\rho _0}{1+d\mathrm{ln}(1+t/T)}.$$
(10)
In this way we can connect directly our model to the experimental data and we can also examine the effect of different dynamical shaking procedures on the compaction process. In this respect, our results are in agreement with those of Nicodemi et al. , who find a similar change of behavior when considering single and multiple tapping. This should clarify to a certain extent the mechanisms of density relaxation in different dynamical situations.

## IV Compaction dynamics with friction

When a frictional force is introduced, as described in Sec. I, some qualitative and quantitative changes are found. In Fig. 12 the time evolution of the density for a fixed value of the friction coefficient $`\mu `$ and several choices of $`\tau `$ is displayed.

A comparison with Fig. 2 shows that for a same value of $`\tau `$ a much smaller asymptotic density is reached, and that a slowing down of the dynamical behavior occurs. As a consequence, larger values of $`\tau `$ are needed to obtain the same final densities as for $`\mu =0`$. In Fig. 13 a comparison of $`\rho (t)`$ between $`\mu =0`$ and $`\mu =0.4`$ for $`\tau =80`$ is illustrated.

It seems that the same functional form Eq. (5) adopted for $`\mu =0`$ can be used to describe the dynamical evolution of the system. In spite of that, for $`\mu 0`$ some slight discrepancies in the tails can be detected for small values of $`\tau `$, not visible at the level of detail of the figure. This means that the reduced probability for horizontal moves affects the compaction mainly in the late stages of the process.

We have also examined the $`\mu `$ dependence of $`\rho (t)`$. The curve for $`\mu =0`$ is far above the others, and furthermore the difference between the curves for low values of $`\mu `$ is more marked. Such differences show up only after an initial transient in which the various curves are practically superimposed one on the others (the curves do not start diverging significantly until the density reaches about $`0.6`$). This is in compliance with what we have just said about the effectiveness of friction in the asymptotic dynamics.

The behavior of the asymptotic density is depicted in Fig. 14 and reveals more interesting aspects. In Fig. 14(a) we have plotted $`\rho _{\mathrm{}}(\tau )`$ for $`\mu =0`$ and $`\mu =0.4`$. We were able to identify a power law regime for small $`\tau `$ also for $`\mu =0.4`$:

$$\rho _{\mathrm{}}(\tau )=(\tau /\tau _l^{(\mu )})^{\delta ^{(\mu )}}$$
(11)
with $`\delta ^{(\mu )}0.062`$ (not very different from $`\delta `$ in Eq. (7)) and $`\tau _l^{(\mu )}6.910^4`$. Thus, saturation density scales in the same way as for $`\mu =0`$, with nearly the same exponent but with a very different relaxation time (more than one order of magnitude larger), i.e. the process is much more sluggish. We did not manage to find the functional behavior for larger $`\tau `$ and it is not very clear whether it is possible to define a crossover shaking time, at least for the values of $`\tau `$ we have considered. For $`\tau >1000`$ $`\rho _{\mathrm{}}`$ increases more steeply than for $`\mu =0`$ and it is likely that one could find another regime for larger $`\tau `$. However computational limitations prevented us to explore in detail the region for $`\tau >10000`$ (since simulations become more and more time consuming as $`\tau `$ increases), and we will address this problem in the future.

Finally in Fig. 14(b) we have analyzed the behavior of $`\rho _{\mathrm{}}`$ vs. $`\mu `$. Interestingly $`\rho _{\mathrm{}}`$ decreases linearly with $`\mu `$ for $`\mu >0.4`$:

$$\rho _{\mathrm{}}(\mu )=AB\mu ,$$
(12)
where $`A0.71`$ and $`B0.073`$. Instead the decrease is more drastic for lower values of $`\mu `$, as noted above. This fact might indicate a change of behavior in the dynamics of the system as a function of the friction coefficient. For this purpose we show in Fig. 15 the stretched exponent $`\beta `$ of Eq. (5) vs. $`\mu `$.

It rapidly decreases for $`\mu <0.4`$, but it is almost constant for larger $`\mu `$. In other terms, the independence of $`\mu `$ of the dynamical evolution reflects itself just in a linear scaling of the final density with $`\mu `$. We hope to elucidate this point more thoroughly in a future work.

## V Conclusions and perspectives

We have proposed and discussed a simplified model for dimer compaction, with the intent to take into account relevant dynamical features, such as horizontal shaking and friction forces between particles, analyzing in detail their effects on the density time evolution. The dynamical behavior of our model is very complex with interesting and novel features. In the absence of friction, for a single tap, the compaction dynamics cannot be interpreted as an inverse logarithmic relaxation but in terms of a stretched exponential law, which is a peculiar characteristic of a continuous range of time constants . Thus, the fact that we do not find a logarithmic law but a stretched exponential relaxation could be a peculiar feature related to the single tap dynamics, since Nicodemi et al.  observe the same kind of behavior in a different model with a similar dynamics. Furthermore, we are not aware of experiments in which a single tap is applied as external perturbation, since they all refer to a tapping sequence. Actually, it turns out that the logarithmic behavior found in experiments is recovered also in our model when a sequential tapping is applied. We therefore might infer that different shaking procedures give rise to intrinsically different compaction behaviors, driven by different dynamical mechanisms. However, it still needs to be clarified whether the logarithmic law is just a good fit to the experimental data or it has a deeper meaning (see for example the discussion in Ref. ).

We also find a crossover in the behavior of the asymptotic density as a function of the shaking time with a slowing down above a critical shaking time value, indicating a possible dynamical transition process. More detailed studies of this point are planned for the future in order to reach a better understanding of the physical mechanism underlying this peculiar behavior.

Similar to what is found in experiments, horizontal shaking favors locally ordered configurations and leads to higher compaction . We believe it would be important to address further this issue by improving our model, allowing rotations of the dimers during the dynamics.

The results obtained by introducing the frictional effect in our model indicate the important role of friction on the compaction dynamics. In particular we find a slowing down of the dynamical behavior which is more evident in the asymptotic regime. The relaxation of the density is still described by a stretched exponential in the presence of friction, with a stretched exponent that becomes constant at large values of the friction coefficient. It would be worthwhile to explore this issue in future works with an improved modeling of the frictional force.