text:
# Optical and Infrared Photometry of the Unusual Type Ia Supernova 2000cx

## 1 Introduction

Astronomers try to understand the universe by looking for patterns in observed phenomena. Often, the patterns themselves are reason to believe in underlying, understandable physical mechanisms, while at other times the exceptions to the rules provide the motivation to expand our conceptions of the physical makeup of cosmic objects. In this paper we present optical and infrared photometry of the very unusual supernova 2000cx. Previous optical data have been presented by Li et al. (2001), who describe SN 2000cx as “unique”.

SN 2000cx was discovered by Yu, Modjaz, & Li (2000) from images taken on 17.5 and 18.4 July 2000 UT as part of the Lick Observatory Supernova Search, using the 0.76-m Katzman Automatic Imaging Telescope (KAIT). This object was located at $`\alpha `$ = 1:24:46.15, $`\delta `$ = +9° 30′ 30$`\stackrel{}{\mathrm{.}}`$9 (equinox 2000.0), which is 23$`\stackrel{}{\mathrm{.}}`$0 west and 109$`\stackrel{}{\mathrm{.}}`$3 south of the nucleus of the S0 galaxy NGC 524. A spectrum taken on 23 July UT with the Nickel 1-m reflector at Lick Observatory (Chornock et al. 2000) revealed that the object was a peculiar Type Ia supernova, resembling SN 1991T a few days before maximum brightness, with prominent Fe III absorption lines near 430 and 490 nm but weak Si II at 612 nm. Optical photometry (Li et al. 2001) revealed that SN 2000cx is different from all known Type Ia SNe and that the light curves cannot be fitted well using the techniques currently available. The pre-maximum rise was relatively fast, similar to SN 1994D, but the post-maximum decline was relatively slow, similar to SN 1991T.

We present optical and infrared photometry of SN 2000cx, initiated at CTIO with the Yale-AURA-Lisbon-Ohio (YALO) 1-m telescope on 19 July 2000 (UT), some 8 days before the time of $`B`$-band maximum. We include data taken with the 0.76-m Manastash Ridge Observatory (MRO) of the University of Washington (also begun on 19 July UT), the CTIO 0.9-m telescope, and the Apache Point Observatory 3.5-m telescope (APO). The calibration of the infrared light curves was primarily made possible with observations made with the Swope 1-m telescope at Las Campanas Observatory.

## 2 Observations

The YALO optical images were obtained using a Loral 2048 $`\times `$ 2048 CCD with a scale of 0$`\stackrel{}{\mathrm{.}}`$30 per pixel, giving roughly a 10′$`\times `$ 10′ field of view. Due to amplifier problems only half of the chip was working before September 2000. On 6 September 2000 this problem was fixed, but it was necessary to change the gain from 6.6 to 3.6 electrons per ADU. The readnoise improved from 14.2 to 11 electrons. The broadband $`JHK`$ imges from the YALO telescope were obtained using a 1024 $`\times `$ 1024 HgCdTe hawaii Array from Rockwell. The filters have over an 80 percent transmittance from 1.171 to 1.322 $`\mu `$m for the $`J`$-band, 1.498 to 1.766 $`\mu `$m for $`H`$ and 2.012 to 2.278 $`\mu `$m for $`K`$. The scale of the infrared CCD is 0$`\stackrel{}{\mathrm{.}}`$20 per pixel, making a total field of view of 3$`\stackrel{}{\mathrm{.}}`$3 $`\times `$ 3$`\stackrel{}{\mathrm{.}}`$3.

The camera used at the CTIO 0.9-m telescope contains the #3 Tektronix 2048 CCD, which is a thinned, anti-reflection coated, back-side illuminated chip with 2K $`\times `$ 2K pixels. The scale and field size are 0$`\stackrel{}{\mathrm{.}}`$40 per pixel and 13$`\stackrel{}{\mathrm{.}}`$5, respectively. The filters used were: $`U`$ (liquid CuSO4), $`\lambda _0`$ = 3575 Å, fwhm = 600 Å; $`B`$, $`\lambda _0`$ = 4202 Å, fwhm = 1050 Å; $`V`$, $`\lambda _0`$ = 5475 Å, fwhm = 1000 Å; $`R`$, $`\lambda _0`$ = 6425 Å, fwhm = 1500 Å; $`I_{KC}`$, $`\lambda _0`$ = 8075 Å, fwhm = 1500 Å. The $`UBVRI_{KC}`$ filters were from the CTIO facility set Tek #1.

Traces of the filters used with the YALO and CTIO 0.9-m telescopes are shown in Stritzinger et al. (2002). However, for our SN 2000cx observations with YALO we did not use the very broad, non-standard $`R`$ filter used by Stritzinger et al. Instead, we used a much narrower, more standard filter (see Appendix).

The APO 3.5-m optical images were obtained using the facility CCD imager SPIcam, which contains a back-side illuminated SITe chip of 2048 $`\times `$ 2048 pixels. 2 $`\times `$ 2 readout was used, giving a scale of 0$`\stackrel{}{\mathrm{.}}`$28 per pixel and a 4$`\stackrel{}{\mathrm{.}}`$78 $`\times `$ 4$`\stackrel{}{\mathrm{.}}`$78 field of view. The APO infrared images were obtained with the 3.5-m telescope using GRIM II, which contains a Rockwell 256 $`\times `$ 256 nicmos3 HgCdTe array. The chip is sensitive from 1 to 2.5 $`\mu `$m with a quantum efficiency of approximately 70 percent, a gain of 4.7 electrons per ADU, and a readnoise of 110 electrons. The $`J`$, $`H`$ and $`K^{}`$ filters transmit at 1.265 $`\pm `$ 0.267, 1.646 $`\pm `$ 0.339, and 2.114 $`\pm `$ 0.343 $`\mu `$m, respectively.

The CCD camera at the MRO 0.76-m telescope uses a Ford Aerospace chip of 1024 $`\times `$ 1024 pixels, with a readnoise of 8 electrons. The scale is 0$`\stackrel{}{\mathrm{.}}`$61 per pixel, giving a 10$`\stackrel{}{\mathrm{.}}`$32 $`\times `$ 10$`\stackrel{}{\mathrm{.}}`$32 field of view. It contains $`UBVRI`$ “Harris” filters.

Instrumental YALO and 0.9-m magnitudes were measured as Point Spread Function (PSF) magnitudes using the daophot package (Stetson 1987, 1990). A transformation equation of the form $`m`$ = f(M,I,X,T), as suggested by Harris, Fitzgerald, & Reed (1981), was used. Here $`m`$ is the observed (i.e. instrumental) magnitude, M is the tabulated magnitude (e.g. from Landolt 1992), I is the tabulated color index, X is the airmass, and T is the time during the night.

Using the observations of the Landolt (1992) standards, we determined zeropoints, color terms and atmospheric terms. This allowed us to calibrate the field stars near SN 2000cx. Photometry of the supernova itself is then tied to the Landolt standards via these field stars. This allows us to derive accurate values for the SN even if it is observed under non-photometric conditions. Once the field star sequences were established, we dropped the airmass term because the differential airmass corrections within a CCD frame are negligible. The time-variable term was dropped because of its demonstrated small contribution to the final photometry (Suntzeff et al. 1999).

Once the MRO and APO images were bias subtracted and flattened, we obtained aperture magnitudes in the iraf environment, using phot within the apphot package, and calibrated the field stars using mknobsfile, fitparams, and evalfit within the photcal package. The transformation equations were configured to produce $`V`$ magnitudes, $`BV`$, $`VR`$, and $`VI`$ colors. The APO infrared data were reduced using iraf scripts written by Alan Diercks, and the IR mosaics were produced using Eugene Magnier’s image processing program mana. We then carried out aperture photometry within iraf.

The mean optical transformation coefficients for the four optical systems described above are to be found in Table 1.

Figure 1 shows a combined $`BVR`$ image taken with the CTIO 0.9-m telescope when SN 2000cx was on the rise. Some nearby field stars are marked. Optical photometry of these stars, obtained from imagery with YALO and MRO, is to be found in Table 2. A comparison of the independent YALO, MRO, and Li et al. (2001) values for the field stars marked as numbers 2, 4, and 8 in Fig. 1 shows good agreement. The range of the mean values for these three stars is 16-29 mmag in $`B`$, 12-17 mmag in $`V`$, 14-30 mmag in $`R`$, and 6-36 mmag in $`I`$, respectively. Thus, photometry of stars can be carried out at the $``$0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$02 mag level.

Absolute calibration in $`J`$, $`H`$, and $`K`$ was done via observations of the SN 2000cx field along with infrared standards of Persson et al. (1998), using the LCO Swope 1-m telescope on five nights in October and November 2001. The mean $`JHK`$ values for field stars 1 and 2 are given in Table 2. Because of the size of the field of view of the IR camera on the APO 3.5-m telescope, only star 1 and SN 2000cx were always on the chip while dithering. To make the fullest use of all the nights when IR data were taken (i.e. photometric and non-photometric nights), we chose to reduce all our IR photometry of SN 2000cx (i.e. from YALO and APO) with respect to field star 1.

Our infrared YALO images were reduced using a package of scripts written by one of us (NBS),222http://www.ctio.noao.edu/ nick/reduction/reduction.html which runs in the iraf environment. This package contain tasks which fill out the file headers with information necessary for subsequent reduction and take care of bias correction, field flattening, masking out bad pixels, and vignetting. Fortunately, field star number 1 and the SN were located out of the zone of vignetting so that we did not have to worry about that at all. A section free of stars was chosen for each night so we could calculate a clean sky level and subtract that from all individual frames.

Tables 3 through 6 contain our optical data from the CTIO 0.9-m, YALO, APO, and MRO telescopes, respectively. In effect, we have four independent optical datasets, with two independent calibrations. The uncertainties of the values in the tables derive from photon statistics and the uncertainties of color corrections, and are to be considered minimum internal errors. A more accurate estimate of the accuracy of our photometry is obtained by fitting fourth order polynomials to the light curves, telescope by telescope and filter by filter, to measure the rms residuals of such fits, under the assumption that the variation of light of the SN is a smooth function for each filter. We have calculated these residuals for our two largest datasets, using data prior to 25 d after the time of $`B`$-band maximum. For the YALO data the internal errors are $`\sigma _B=\pm `$ 14; $`\sigma _V=\pm `$ 34; $`\sigma _R=\pm `$ 41; and $`\sigma _I=\pm `$ 71 mmag. For the MRO data the internal errors are $`\sigma _B=\pm `$ 6; $`\sigma _V=\pm `$ 20; $`\sigma _R=\pm `$ 43; and $`\sigma _I=\pm `$ 51 mmag.

Figure 2 shows our optical data. Fig. 3 is the same as Fig. 2, but with the addition of the KAIT and Wise Observatory data given by Li et al. (2001), plus data obtained by Jha (2002) with the 1.2-m telescope at the Fred L. Whipple Observatory at Mt. Hopkins, Arizona.

There are systematic differences between different subsets of the data. With the sense of $`\mathrm{\Delta }`$ “YALO minus MRO”, at the time of $`B`$-band maximum $`\mathrm{\Delta }B=`$0.04, $`\mathrm{\Delta }V=`$0.05, $`\mathrm{\Delta }R=`$0.04, $`\mathrm{\Delta }I`$ = +0.05 mag. At $`t`$ = 17.4 d, from Tables 4 and 6 we find $`\mathrm{\Delta }B=`$0.02, $`\mathrm{\Delta }V=`$0.08, $`\mathrm{\Delta }R=`$0.15, $`\mathrm{\Delta }I`$ = +0.16 mag. The disagreement in $`V`$, $`R`$, and $`I`$ is large.

Now, with the sense of $`\mathrm{\Delta }`$ “CTIO 0.9-m minus MRO”, for the data within 7 days of maximum we find $`\mathrm{\Delta }B=`$0.02, $`\mathrm{\Delta }V=`$0.04, $`\mathrm{\Delta }R=`$0.05, $`\mathrm{\Delta }I`$ = 0.00 mag. This is not a significant improvement on the YALO vs. MRO situation at maximum. However, YALO photometry vs. CTIO 0.9-m data within a week of maximum is in agreement at the 0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$02 level or better for $`B`$, $`V`$, and $`R`$. In spite of systematic differences between datasets obtained with different telescopes and different filters, at least we know the apparent magnitudes at maximum to $`\pm `$ 0.03 mag.

The systematic differences in the photometry are undoubtedly due to differences in the actual filters used, coupled with the non-stellar spectral energy distribution of the SN. Very late-time $`I`$-band data are particularly discordant, more than 0.7 mag (KAIT vs. FLWO). Though some of the optical spectra discussed by Li et al. (2001) cover the full $`I`$-band (at 2, 6, 7, 32, 42 days after T($`B_{max}`$) and later), to reconcile all data of SN 2000cx taken with different telescopes is beyond the scope of this paper.333Stritzinger et al. (2002) attempted to correct data of SN 1999ee obtained with multiple telescopes by calculating “S-corrections” using the spectra of this object plus knowledge of the filter transmission curves, quantum efficiencies of the chips, and appropriate atmospheric transmission functions. In the end they decided to leave their data uncorrected. However, more recently Krisciunas et al. (2003) derived corrections to their photometry of SN 2001el and applied the filter corrections to the $`BVJHK`$ data. This in particular solved problems with the $`B`$-, $`V`$-, $`J`$-, and $`H`$-band data. However, they found that $`R`$-band corrections were not really necessary, and that applying $`I`$-band corrections actually made datasets obtained with different telescopes much more discordant.

Tables 7 and 8 contain the infrared data from YALO and APO, respectively. Table 9 has corrections to place the APO data (which are fewer in number) on the YALO filter system, using the method of Stritzinger et al. (2002) and Krisciunas et al. (2003). To our knowledge, only two infrared spectra of SN 2000cx itself exist (Rudy et al. 2002), and these were taken 6 and 7 days, respectively, before the observed date of $`B`$-band maximum. To calculate the filter corrections for the APO data of SN 2000cx we used the infrared spectra of SN 1999ee (Hamuy et al. 2002). In Figs. 4 and 5 we show our infrared photometry without, and with, the filter corrections. The agreement of the $`J`$-band datasets is clearly better with the corrections.

Finally, in Table 10 we give the times of maximum, as observed by us in the different bands. Our time of $`B`$-band maximum is 0.3 d later than that found by Li et al. (2001), well within the uncertainty. While they found that the $`V`$-band maximum occurred 2.1 d after $`B`$-band maximum, we find a lag of only 1.2 d. However, this difference is not statistically significant.

Jha (2002) also obtained $`UBVRI`$ photometry of SN 2000cx, though with a large gap in time after maximum. Late time photometry (10 July 2001) with $`HST`$, using the $`F675W`$ and $`F814W`$ filters, is given by Li et al. (2002). Altogether, the photometric database of SN 2000cx amounts to 642 data points. To our knowledge it is the largest dataset ever obtained for a Type Ia SN.

## 3 General Discussion

A vast majority of the light curves of Type Ia SNe can be fit into a classification scheme in which the shape of the light curves is correlated with the intrinsic luminosity of the SN (Phillips 1993, Riess, Press, & Kirshner 1996, Perlmutter et al. 1997, Phillips et al. 1999). Li et al. (2001) remark that the Multi-color Light Curve Shape (MLCS) fit to their data for SN 2000cx is the worst fit they have ever seen. SN 2000cx was a reasonably fast riser, but a slow decliner. That would mean that the “stretch factor” used with the method of Perlmutter et al. (1997) would give a different value prior to maximum compared to after maximum. Since the stretch factor in $`B`$ and $`V`$ is related to the intrinsic luminosity, how does one determine the intrinsic luminosity of this SN?

NGC 524, the host galaxy of SN 2000cx, has had its distance modulus measured via the Surface Brightness Fluctuation (SBF) method. Tonry et al. (2001) obtain $`mM`$ = 31.90 $`\pm `$ 0.20 for the distance modulus. Li et al. (2001) obtained $`mM`$ = 32.53 $`\pm `$ 0.35 using MLCS. Using the corrected recession velocity of NGC 524 in the Local Group frame of 2192 km s-1 and a Hubble constant of 74 km s-1 Mpc-1, in agreement with the $`HST`$ Key Project value (Freedman et al. 2001), we obtain a distance modulus of 32.36 mag.

SN 2000cx was located in the outer regions of an early-type galaxy. Since early-type galaxies contain minimal amounts of dust, and the $`BV`$ colors of SN 2000cx were particularly blue, we explicitly assume that this SN was unreddened in its host. We shall correct for reddening due to dust in our Galaxy (see below).

Krisciunas et al. (2003) showed that the $`H`$-band absolute magnitudes of Type Ia SNe 10 days after the time of $`B`$-band maximum appear to be a flat function of the decline rate $`\mathrm{\Delta }`$m15($`B`$). For a sample of 9 objects with $`\mathrm{\Delta }`$m15($`B`$) $`<`$ 1.3 they found a mean $`H`$-band absolute magnitude of $``$17.91 $`\pm `$ 0.05. The observed $`H`$-band magnitude of SN 2000cx 10 days after $`B`$-band maximum is 14.64 $`\pm `$ 0.04. From the Galactic reddening maps of Schlegel, Finkbeiner, & Davis (1998), we estimate that the color excess E($`BV`$) = 0.082 mag in the direction of SN 2000cx. Using the interstellar extinction model of Rieke & Lebofsky (1985), we estimate that the $`H`$-band extinction, due only to the effect of dust in our Galaxy, is 0.04 mag. Assuming that the absolute magnitude of SN 2000cx is also $``$17.91, given its corrected $`H`$-band magnitude of 14.60, we obtain a distance modulus $`mM`$ = 32.51. Thus, the distance moduli from MLCS, Hubble’s Law, and the $`H`$-band analysis are in reasonable agreement, with an unweighted mean value of 32.47 $`\pm `$ 0.05, corresponding to a distance of 31 $`\pm `$ 1 Mpc.

Given the observed $`V`$-band maximum of 13.25, A$`{}_{V}{}^{}`$ 3.1 $`\times `$ 0.082 = 0.25, and a distance modulus of 32.47, the implied $`V`$-band absolute magnitude of SN 2000cx is $``$19.47, which is comparable to the mean of MV of the slow decliners (Krisciunas et al. 2003, Fig. 13).

However, Ajhar et al. (2001) have shown that there is excellent agreement between distance determinations using SBF and other methods. The host of SN 2000cx is actually a bit close ($`cz<3000`$ km s-1) to derive its distance via Hubble’s Law. The MLCS distance can also be doubted because the light curves cannot be fit well using MLCS templates. Also, on the basis of a larger sample of objects, it may turn out that $`H`$-band absolute magnitudes do show some kind of decline rate relation. In that case, should we use the decline rate of SN 2000cx or some modified value that takes into account its different rise and decline rates?

Let us assume that the Tonry et al. (2001) distance modulus of $`mM`$ = 31.90 $`\pm `$ 0.20 is correct. On an $`H_0`$ = 74 scale this would be $`mM`$ = 31.84 mag. Adopting AB = 0.34, AV = 0.25, and AI = 0.12 mag and the maximum magnitude values given in Table 10, we obtain absolute magnitudes of MB = $``$18.76, MV = $``$18.84, MI = $``$18.31, with uncertainties of $`\pm `$ 0.20 mag. These correspond to $`\mathrm{\Delta }`$m15($`B`$) in the range 1.4-1.7 (Krisciunas et al. 2003, Fig. 13). MH ($`t`$ = 10 d) would be $``$17.24, comparable to that of SN 2000bk, which was a fast decliner with $`\mathrm{\Delta }`$m15($`B`$) = 1.63.

Krisciunas et al. (2001, Figs. 16 and 17) devised a quantitative measure of the strength of the $`I`$-band secondary hump common to Type Ia SNe, namely the mean flux-with-repect-to-maximum from 20 to 40 d after the time of $`B`$-band maximum. They found that 90 percent of Type Ia SNe have values of $``$ I $`_{2040}`$ that are well correlated with the $`B`$-band decline rate $`\mathrm{\Delta }`$m15($`B`$). There were two exceptions to the rule, showing that there can be objects with identical decline rates in $`B`$ and $`V`$ but much stronger or weaker secondary humps in $`I`$. For SN 2000cx we find that $``$ I $`_{2040}`$ = 0.35. This would imply $`\mathrm{\Delta }`$m15($`B`$) $``$ 1.7 like the fast decliners SNe 1992bo or 1993H. But SN 2000cx is a slow decliner, with $`\mathrm{\Delta }`$m15($`B`$) = 0.93. If we were to add a point in Fig. 17 of Krisciunas et al. (2001) corresponding to SN 2000cx, it would be the most discrepant point in the graph.

Li et al. (2001) point out that the stretch factor for the pre-maximum data points ($`t=8`$ to 1 day) corresponds to $`\mathrm{\Delta }`$m15($`B`$) = 1.64 $`\pm `$ 0.02. In various ways SN 2000cx masquerades as a fast Type Ia SN and also as a slow one.

One of the patterns exhibited by many unreddened Type Ia SNe is that their $`BV`$ colors from 30 to 90 d after the time of $`V`$-band maximum follow a particular linear trend, whether they are fast-decliners or slow-decliners (Lira 1995, Phillips et al. 1999). Li et al. (2001) and Cuadra et al. (2001) noted that SN 2000cx was roughly 0.2 mag bluer than the Lira line. In Fig. 6 we show the $`BV`$ data from YALO, the CTIO 0.9-m, APO, and MRO telescopes. The data have been dereddened by 0.082 mag to account for dust in our Galaxy (Schlegel et al. 1998). The data after JD 2,451,784 are, on average, 0.217 mag below the Lira line. We assume, because SN 2000cx occurred in the outer regions of an early-type galaxy, that it suffered no host reddening. But if its light was affected by dust in the host galaxy, then the points in Fig. 6 should be displaced even further below the Lira line. Truly, the $`BV`$ colors of this SN are unusual.

In Figs. 7, 8, and 9 we show the $`JHK`$ light curves of SN 2000cx along with data of other objects: SNe 1999aw (Strolger et al. 2002), 2001el (Krisciunas et al. 2003), 1999ac (Phillips et al. 2002, 2003), 2000bk (Krisciunas et al. 2001), and 1986G (Frogel et al. 1987). The light curves are ordered from top to bottom by the decline rate $`\mathrm{\Delta }`$m15($`B`$).

We the note the extremely deep dip in the $`J`$-band in the 20 days after maximum. We also note the similarity of the SN 2000cx $`J`$-band data with that of SN 2000bk, a fast decliner. The latter had an earlier secondary $`J`$-band maximum, however. We also note, whereas most Type Ia SNe have relatively flat $`H`$\- and $`K`$-band light curves in the 20 days after maximum, SN 2000cx showed decreasing flux at this epoch.

Krisciunas et al. (2000, 2001, 2003) found that Type Ia SNe which are mid-range decliners appear to exhibit uniform $`Vminus`$ infrared color curves from about one week before maximum until three or more weeks after maximum. They asserted that the unreddened loci can be used to derive the total extinction suffered by the light of a Type Ia SN on its way to our viewing location in the Galaxy. They found that $`VH`$ and $`VK`$ were particularly well behaved, and modelling by Höflich given in Krisciunas et al. (2003) confirms this observational result from a theoretical standpoint. Krisciunas et al. (2000, 2001) also noted that fast decliners and slow decliners have different unreddened loci. The slow decliners have bluer loci, and the fast decliners have redder loci.

Using data of the slow decliners SNe 1999aa (Krisciunas et al. 2000), 1999aw (Strolger et al. 2002), and 1999gp (Krisciunas et al. 2001), which appear to be unreddened in their hosts, we can correct the $`VJ`$, $`VH`$, and $`VK`$ colors for the effect of dust in our Galaxy (Schlegel et al. 1998) and construct unreddened loci for slowly declining Type Ia SNe. We also used the optical data of SN 1999ee (Stritzinger et al. 2002) along with some unpublished IR data of SNe 1999ee and 2001ba taken at Las Campanas Observatory and CTIO to constrain the shape of the unreddened loci.

Following Krisciunas et al. (2000), let $`t`$ equal the number of days since the time of $`B`$-band maximum, and Tc be some “crossover time” when the slope in the color curve has a sudden change. We consider only the data from $`9<t<27`$ d and construct the simplest form of unreddened loci that can fit the data. For $`VJ`$ we have:

$$VJ=a_0+a_i+b_1(tT_c)+c(tT_c)^2;(t<T_c).$$
(1)
$$VJ=a_0+a_i+b_2(tT_c);(tT_c).$$
(2)
Each supernova has its own ai value, but the objects presumed to be unreddened and which are used to establish the unreddened locus will have $``$ a$`{}_{i}{}^{}`$ = 0. Reddened objects have a$`{}_{i}{}^{}>0`$. Under the assumption that the color curves of a group of unreddened supernovae are validly parameterized by a uniform locus, and that the reddened objects exhibit the same locus simply shifted to the red, the goal is to use all the data for a given color index to solve for the crossover time, the zeropoint, the slopes, the second order coefficient, and the color excesses of each reddened object.

In Fig. 10 we show the unreddened $`VJ`$ color locus for slowly declining Type Ia SNe. Using the $`BV`$ color excesses given by Schlegel et al. (1998) and E($`VJ`$) = 2.223 E($`BV`$) from a standard model of Galactic dust (Rieke & Lebofsky 1985), we have subtracted $`VJ`$ color excesses of 0.089, 0.072, and 0.124 mag, respectively, from the observed $`VJ`$ colors of SNe 1999aa, 1999aw, and 1999gp. We assume that these three objects were unreddened in their host galaxies. For SN 1999ee we derive E($`VJ`$) = 0.658 $`\pm `$ 0.022. This has been subtracted from the SN 1999ee data shown in Fig. 10. The $`V`$ near IR color loci given by Krisciunas et al. (2000) were not extrapolated beyond $`t`$ = 27 d, and we have truncated our $`VJ`$ locus at that epoch, though on the basis of SN 1999aw we might have extended the right-hand line further.

The $`VJ`$ photometry of slowly declining Type Ia SNe gives a crossover time Tc = 14.43 $`\pm `$ 0.33 d. We find $`a_0`$ = $`1.96`$, $`b_1`$ = $``$0.07480 $`\pm `$ 0.00747 mag d-1, $`c`$ = $``$0.00015 $`\pm `$ 0.00030 mag d-2, and $`b_2`$ = 0.06867 $`\pm `$ 0.00500 mag d-1. Whereas the $`c`$ term was necessary for the mid-range decliners studied by Krisciunas et al. (2000), we found that it was not statistically significantly different than zero for the slow decliners. We list it here only for reasons of completeness.

Given the optical color excess E($`BV`$) = 0.082 for SN 2000cx and the interstellar extinction law quantified by Rieke & Lebofsky (1985), we estimate the following color excesses: E($`VJ`$) = 0.182; E($`VH`$) = 0.210; and E($`VK`$) = 0.226.

In Fig. 11 we show the dereddened $`VJ`$ colors of SN 2000cx, along with the unreddened color locus described above (the solid line). We also show for reference (as a dashed line) the unreddened color locus for mid-range decliners given by Krisciunas et al. (2000). The big dip observed in the $`J`$-band leads to a very blue $`VJ`$ color for SN 2000cx at $`t`$ = 10 d. However, the $`VJ`$ colors match the unreddened locus for slow decliners overlapping the time of maximum light. We note that conversion of the $`J`$-band data to the filter system of Persson et al. (1998) would make our data up to 0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$07 fainter at maximum light, assuming that the IR spectral evolution of SN 2000cx was similar to SN 1999ee (see Krisciunas et al. 2003, Fig. 6). This would make the $`VJ`$ data shown on the left hand side of Fig. 11 up to 0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$07 bluer.

A similar analysis of the $`VH`$ colors of SNe 1999aw, 1999gp, and 1999ee gives a crossover time Tc = 5.71 $`\pm `$ 0.27 d, $`a_0`$ = $`1.60`$, $`b_1`$ = $``$0.06872 $`\pm `$ 0.00304 mag d-1, $`c`$ 0, and $`b_2`$ = 0.07502 $`\pm `$ 0.00262 mag d-1. In this case we have subtracted E($`VH`$) = 0.083 mag from the observed colors of SN 1999aw and 0.143 mag from the data of SN 1999gp, to account for the effect of dust in our Galaxy. We derived a total color excess of E($`VH`$) = 0.795 $`\pm `$ 0.023 for the reddening of SN 1999ee. The E($`VJ`$) and E($`VH`$) color excesses of SN 1999ee are consistent with one unique value (rather than two disjoint values) of AV = 0.944 $`\pm `$ 0.061, which is a good consistency check relating to the assumption that there exist uniform color loci.

In Fig. 12 we show a corresponding plot for dereddened $`VH`$ colors of SN 2000cx, with the unreddened loci for slow decliners (solid line) and mid-range decliners (dashed line) shown. For $`VH`$ there is good agreement between the SN 2000cx data and the unreddened locus for slow decliners from $`t`$ = 10 to 27 d. Conversion of the $`H`$-band data to the filter system of Persson et al. (1998) would change the data from $`10<t<27`$ d by up to 0$`\stackrel{\mathrm{m}}{\mathrm{.}}`$06, making them brighter, and the $`VH`$ colors redder. Once again, this assumes that the IR spectral evolution of SN 2000cx was similar to that of SN 1999ee.

In Fig. 13 we show the $`VK`$ data of SN 2000cx, corrected only for the reddening due to our Galaxy, along with dereddened data of SNe 1999aa, 1999ee, 1999gp, and 2001ba. The $`VK`$ colors of SN 2000cx are quite similar to these other slow decliners, both pre- and post-maximum. From $`10<t<21`$ d SN 2000cx is the bluest object in this color index.

## 4 The Bolometric Behavior of SN 2000cx

The wide wavelength coverage of the $`UBVRIJHK`$ broadband magnitudes allows us to construct ultraviolet/optical/near-infrared “uvoir” bolometric light curves. Only a few papers have been published trying to estimate the bolometric light curves of Type Ia supernovae. (See Leibundgut & Suntzeff 2003 for a summary.) The calculated luminosities are not the true bolometric luminosities, but represent the fraction of the gamma rays produced in the radioactive decays of the synthesized nuclides that are thermalized in the expanding debris nebula. As shown by Leibundgut & Pinto (1992) and Leibundgut (2000), a significant fraction of the gamma rays leak out of the supernova debris going from 10 percent at the time of $`B_{max}`$ to over 50 percent 40 days after maximum.

Suntzeff (1996) used ultraviolet spectra, optical $`UBVRI`$, and near-infrared $`JHK`$ data to estimate accurate bolometric fluxes. Because of the limited infrared data at the time, only a few points on the bolometric light curves could be accurately calculated. Vacca & Leibundgut (1996) and Contardo, Leibundgut, & Vacca (2000) obtained more complete uvoir bolometric light curves by integrating optical broad-band magnitudes. Applying Arnett’s law (Arnett 1982) to the peak bolometric luminosities, Vacca & Leibundgut (1996) and Contardo et al. (2000) found a range of more than a factor in 10 in the 56Ni masses for a group of nearby Type Ia SNe. Cappellaro et al. (1997) used a $`V`$ magnitude with a bolometric correction to study the gamma-ray trapping in the late-time light curves, which also showed a significant range in 56Ni masses.

We have integrated the broad-band magnitudes for SN 2000cx using a simple trapezoidal rule and a conversion of broad-band to monochromatic fluxes from Suntzeff & Bouchet (1990). We have added extrapolations to the ultraviolet from the $`U`$ filter and to the mid-infrared from the $`K`$ or $`H`$ filter using an extrapolation scheme discussed by Suntzeff (2003). These extrapolations only add $`2`$ percent for the missing infrared flux, and less than 10 percent for the ultraviolet after maximum light. In Fig. 14 we plot the uvoir bolometric luminosity for SN 2000cx. For comparison, we also plot similar data for SN 2001el and SN 1999ee taken from Suntzeff (2003). (Note that the data have been offset slightly for plotting purposes.) The SN 1999ee and SN 2001el data were taken from Stritzinger et al. (2002) and Krisciunas et al. (2003) and are meant to represent “normal” Type Ia supernovae with $`\mathrm{\Delta }`$m15($`B`$) of 0.94 and 1.13, which are similar to the $`\mathrm{\Delta }`$m15 value of SN 2000cx. We have assumed the following distance moduli and $`BV`$ reddening, which are based on a distance scale of $`H_0=74`$ km s-1 Mpc-1: SN 2000cx, (32.47,0.082); SN 1999ee, (33.21,0.30); SN 2001el, (31.26, 0.21).444A systematic error of $`\mathrm{\Delta }`$m in the distance modulus corresponds to a change in $`log`$(L) of 0.4 $`\times \mathrm{\Delta }`$m. If the true distance modulus of SN 2000cx is 31.90, then we must shift the bolometric light curve of SN 2000cx in Figs. 14 and 15 down by 0.228 mag. The bolometric light curves have been plotted relative to the time of peak bolometric luminosity (not $`B_{max}`$) which we estimated from the curves.

The only remarkable difference between SN 2000cx and the normal Type Ia’s as seen in Fig. 14 is the lack of the flux excess around day 30. This flux excess, which is associated with the secondary maxima in $`IJHK`$ and a “shoulder” in $`R`$ was noted by Suntzeff (1996) and Contardo et al. (2000). Since the energy source for the bolometric luminosity and the optical depth to gamma rays are monotonically declining at that date, such an inflection in the bolometric luminosity must be associated with cooling mechanism and not the energy input to the nebula.

In Fig. 15 we can see the differences in the bolometric light curves more clearly. Here we plot smoothed representations of the bolometric light curves, and at the bottom of the panel, the bolometric light curve of SN 2000cx with respect to SN 1999ee and SN 2001el. It can be seen that SN 2000cx rises to maximum more quickly and falls more quickly (within 10 days of maximum), but the effect is rather small. However, starting around day 25, SN 2000cx suddenly declines rapidly compared to these two SNe, reaches a maximum flux deficit around day 35, and then increases in brightness slightly to day 50. Fig. 15 gives the impression that the post-maximum flux enhancement of SN 2000cx was weaker and earlier than the comparison Type Ia events.

Two explanations for the secondary maximum have been published. One explanation for the secondary $`IJHK`$ maxima has been given by Pinto & Eastman (2000a, 2000b) who note that this flux enhancement is due to a rapid change in the flux mean opacity. After maximum light, the thermalized energy input to the light curve from the radioactive nuclides is less than the observed luminosity, implying qualitatively that the post-maximum luminosity is powered by a reservoir of previously trapped radiation. If the post-maximum opacity decreases due to a drop in the effective temperature, the diffusion times drop and the trapped energy escapes more rapidly leading to a pause in the rapid luminosity decline. This bolometric flux excess appears in the redder colors because the opacities are very low and there are ample emission sources such as Fe II and Ca II.

A similar explanation has been given by Höflich (1995) and Höflich, Khokhlov, & Wheeler (1995). They note that the infrared luminosity can be roughly approximated by the Rayleigh-Jeans limit, namely:

$$(L_{IRcolor})R_{ph}^2T_{eff}D_{IR}$$
(3)
Here $`R_{ph}`$, $`T_{eff}`$, and $`D_{IR}`$ are the photospheric radius, the effective temperature, and a dilution factor appropriate for a scattering dominated atmosphere. They argue that $`D_{IR}`$ is a slowly varying function during this epoch. They show that after maximum, $`T_{eff}`$ drops as the energy source switches from 56Ni to 56Co which starts the steep post-maximum decline. In many of their models however, the photospheric radius $`R_{ph}`$ increases as it is dragged out by the expansion of the debris. Depending on the rate of expansion of the photosphere, this can cause the product in the equation for $`L_{IRcolor}`$ to increase. The appearance of an expanding photosphere can only be maintained if the opacities stay high. In these conditions, as the photosphere is pulled outward with the debris, the luminosity will increase if $`T_{eff}`$ does not drop dramatically. This high rate of cooling, which greatly exceeds the energy input from the radioactive nuclides, cannot be maintained indefinitely, and at some point the opacities will begin to drop so rapidly that apparent photospheric radius will also being to recede quickly causing a sudden decreace in luminosity.

According to the models discussed above, this would support the hypothesis that SN 2000cx is a sub-luminous event. The sub-luminous SNe 1992A, 1992bo, and 1991bg (Contardo et al. 2000; see their Figure 5) also have weak or no secondary flux enhancements and SN 2000cx clearly belongs to this class.

Unfortunately, given the large number of unknown parameters in the models for the explosions of Type Ia SNe, this is not a strong conclusion. Pinto & Eastman (2000b) note that the secondary maximum is a sensitive function of how the much of the radioactive nuclides are mixed in the expanding debris from the core. An unusual mixing event bringing relatively more 56Ni out from the core qualitatively could account for the lack of the secondary flux enhancement independent of the intrinsic luminosity.

Returning to the empirical bolometric light curve, we find two other aspects of the morphology of the bolometric light curve which point to this event being sub-luminous. Contardo et al. (2000) showed that the decline rate between days 50 and 80 for their sample of Type Ia supernovae was $`2.6\pm 0.1`$ mag per 100 days, except for most sub-luminous event in their sample, SN 1991bg, which had a decline rate of 3.0 mag per 100 days. SN 2000cx declined at 2.9 mag per 100 days during this time period. If one looks at their Fig. 5, one can also see that SN 1991bg also stands out in the peak-to-tail luminosity difference. In Fig. 16 we plot the difference in luminosity between the peak luminosity and 90 days after peak for the sample studied by Contardo et al. (2000). Evidently, this luminosity difference is correlated with $`\mathrm{\Delta }`$m15($`B`$), with the fainter SNe having larger luminosity differences. The observed value of this luminosity difference in SN 2000cx of 1.63 associates it with the sub-luminous group. This is also seen in Fig. 15 where this SN is compared to SNe 2001el and 1999ee. In that figure, the luminosity on the exponential tail for days beyond 50 are underluminous by about 0.2dex with respect to SN 2001el, which has a value of $`\mathrm{\Delta }`$m15($`B`$) of 1.13.

By day 50 or so, the energy deposition in a typical Type Ia supernova due to the thermalization of gamma rays occurs in regions which are optically thin in the optical and near-infrared (Pinto & Eastman 2000a). Thus, the luminosity at this epoch responds rapidly to the input energy source, which at this time is 56Co. At maximum light, according to Arnett’s law, the bolometric luminosity equals the instantaneous energy input from the radioactive nuclides. The ratio of these two luminosities should then be independent, to first order, of the mass of 56Ni synthesized.

The larger luminosity difference between the peak and 90 days after peak is thus caused by a smaller optical depth to gamma rays at late times. This could be caused, for instance, by positron escape as discussed by Milne, The, & Leising (2001) but the modelling shows minimal effects of positron trapping at this epoch due to the short positron lifetimes which approximate in situ deposition of positron energy. A lower optical depth to gamma rays could also be caused by an asymmetric mass distribution of the ejecta. Perhaps the simplest way to lower the optical depth is to increase the kinetic energy due to the explosion.

It is not unreasonable that the initial kinetic energy may be only vaguely related to the mass of 56Ni synthesized. Pinto & Eastman (2000a; see their Figure 4), showed that a larger kinetic energy will lead to a more rapid decline in the mass column depth, and produce a narrower bolometric light curve and increased the peak-to-tail luminosity difference. In a subsequent article, Pinto & Eastman (2001) note that the explosion kinetic energy is not necessarily a function of the 56Ni mass, since the total energy of the burning of a C/O mixture is nearly the same if it burns to the Si group or to 56Ni. The thermal energy liberated at the explosion does go entirely to the kinetic energy, but that kinetic energy may not be strongly coupled to the 56Ni, which powers the subsequent light curve. Li et al. (2001) found that SN 2000cx did have very high sulfur and silicon velocities compared to SN 1994D, and also concluded that SN 2000cx may have had a larger than typical kinetic energy.

Thus, if the SBF distance is correct, a working hypothesis to explain the bolometric behavior of SN 2000cx is that this supernova is an underluminous event with higher than normal kinetic energy. We should not be too forceful in stressing this conclusion however. This event was clearly unusual in its color evolution, and simple morphological arguments may fail if model parameters, such as the extent of mixing of the radioactive nuclides or the symmetry of the explosion, are uncoupled from the other fundamental parameters such as the amount of 56Ni synthesized. None of the evidence here clearly points to a shorter or longer distance to this supernova. Obviously, we need a better distance to NGC 524 to resolve this question.

## 5 Conclusions

SN 2000cx, the brightest supernova discovered in the year 2000, occurred in the unobscured outer regions of an early-type galaxy and was well observed with multiple telescopes, allowing us to compile a datset of unprecedented size. While many Type Ia SNe have light curves that follow patterns that are now well established, SN 2000cx did not conform to these patterns.

SN 2000cx was a reasonably fast riser in $`B`$ and $`V`$. From the pre-maximum photometry Li et al. (2001) obtain a stretch factor that corresponds to $`\mathrm{\Delta }`$m15($`B`$) = 1.64 $`\pm `$ 0.02. Based solely on its weak $`I`$-band secondary hump, we would have predicted $`\mathrm{\Delta }`$m15($`B`$) $``$ 1.7. If the distance modulus based on Surface Brightness Fluctuations of the host galaxy is correct (Tonry et al. 2001), the corresponding absolute magnitudes in $`BVIH`$ are comparable to fast decliners, with $`\mathrm{\Delta }`$m15($`B`$) in the range 1.4 to 1.7.

However, SN 2000cx was a slow decliner, with $`\mathrm{\Delta }`$m15($`B`$) = 0.93. Its pre-maximum spectrum showed strong Fe III and weak Si II, like other slow decliners such as SN 1991T. Its $`VK`$ color evolution, both pre- and post-maximum, was very similar to that of other slow decliners.

The bolometric behavior of SN 2000cx, when compared to the normal SNe 1999ee and 2001el, showed that this SN rose and fell from maximum light more rapidly, and that the magnitude difference between peak brightness and 90 days past peak was larger than normal. This behavior is consistent with the higher kinetic energies seen at maximum light and reported by Li et al. (2002), but it can also be explained by this event being sub-luminous.

The distance modulus of SN 2000cx is somewhat problematic. We note, however, that MLCS (Riess et al. 1996, 1998), Hubble Law’s (with H0 = 74 km s-1 Mpc-1), and the method that uses the $`H`$-band absolute magnitude at $`t`$ = 10 d (Krisciunas et al. 2003) give just about the same value. The SBF method gives a distance modulus roughly 0.6 mag smaller.

Given: 1) the accuracy of the photometry of SN 2000cx at maximum ($`\pm `$ 0.03 mag); 2) the host extinction of SN 2000cx is minimal (or zero); 3) the uncertainty of the Galactic extinction correction is also small; and 4) the light curves cannot be fit by templates based on other objects, confidently placing this SN in a Hubble diagram depends significantly on a direct measure of the host galaxy’s distance (such as with the SBF method) or on an absolute magnitude derived from an explosion model that can match the many unusual observed facts.

Support for Proposal Number GO-07505.02A, GO-08177.6, and GO-08641.07A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This paper is based, in part, on observations obtained with the Apache Point Obsevatory 3.5m telescope, which is owned and operated by the Astrophysical Research Consortium. We thank W. D. Li for sharing the KAIT data considerably ahead of publication. We thank M. R. Garcia, L. Clark, D. Hoard, J. Alfonso, and K. Vivas for obtaining some of the data at CTIO. C. Stubbs, G. Miknaitis, and E. Bergeron helped with data acquisition for some APO observations. D. Edgeworth helped with some of the MRO observations. The APO infrared data were reduced using software written in part by A. Diercks and E. Magnier.

## Appendix A The $`R`$-band filter used at YALO

In Fig. 17 we show the filter transmission function of the $`R`$-band filter used with the YALO 1-m telescope for the observations presented here. This is a much more standard filter than that used by Stritzinger et al. (2002) for their observations of SN 1999ee using the same telescope and camera. Also shown in the diagram is the effective throughput in $`R`$, made up of the filter transmission function multiplied by an atmospheric transmission function, the quantum efficiency of the chip, and two aluminum reflections. Both functions shown in Fig. 17 include a 280 Å shift to the blue, which comprises the 310 Å shift suggested by the manufacturer, taking into account the use of the filter when cooled, and a 30 Å shift back to red. This smaller shift was necessary to match the color term obtained for $`R`$-band photometry based on synthetic photometry with the actual color term derived from the observations of Landolt (1992) standards.


format: text

text:
# A phenomenology for multiple phases in the heavy fermion skutterudite superconductor PrOs4Sb12

## I Introduction

Recently, heavy fermion superconductivity was discovered at 1.85 $`K`$ in the skutterudite compound PrOs4Sb12 with a cubic (tetrahedral) crystal symmetry Th. Intriguingly, it was indicated that there are multiple superconducting phases in this material, like in the Uranium compounds. Two specific heat jumps were observed at $`T_{c1}=1.85K`$ and $`T_{c2}=1.75K`$ in the absence of the magnetic field. Moreover, the thermal conductivity measurement in a magnetic field rotated relative to the crystal axes reveals a novel change in the symmetry of the gap function in the $`ab`$-plane. The result was that the gap function has four point-node-like structure (four dips) in the  and  directions on the Fermi surface in the high field and high temperature region (A-phase), and two dips in the  direction in the low field and low temperature region (B-phase). The proposed phase diagram is given in Fig. I. Unfortunately, the gap structure along the $`c`$-axis remains unresolved within in the presently available data, and one of the purpose of this paper is to consider this problem. The presence of the point node like structures appears to be consistent with the power law behavior of the specific heat and the nuclear-spin-lattice relaxation at low temperatures (See also Ref. ). Another evidence for the double superconducting transition is the thermal expansion measurement along  direction; two clear jumps are observed at around $`T_{c1}=1.85K`$ and $`T_{c2}=1.75K`$ in the absence of the magnetic field (Fig. I).

In this paper, the symmetry of the gap functions for A- and B-phase are discussed by use of a phenomenological Ginzburg-Landau (GL) theory based on the tetrahedral crystal symmetry Th. At present, there are no information as to whether the superconductivity is spin singlet or spin triplet pairing. Only the singlet pairing case will be discussed in this paper. The gap structure along the $`c`$-axis is considered. We find that the gap function appears to have two additional dips along the $`c`$-axis in A-phase (in total there are six dips, along the ,  and  directions) , while no additional dips are inferred for the B-phase resulting in a toal of two dips, along the  direction) . The A-phase gap is essentially an ”anisotropic $`s`$-wave” state and the B-phase gap is an ”anisotropic $`s+id`$-wave” state. The B-phase gap breaks the cubic crystal symmetry and time reversal symmetry spontaneously, from which one can expect a crystal deformation with cubic symmetry breaking to occur in this phase. To verify this, the phenomenological strain-order parameter coupling is taken into account and the thermal expansions and isothermal elastic constants are calculated. One finds that these values become anisotropic in B-phase.

The pairing mechanism for this superconductivity should be also interesting. It is widely believed that the mechanism for the heavy fermion superconductivity in the Uranium compounds comes from magnetic interactions, but PrOs4Sb12 is a non-magnetic compound. Several experiments indicate that the ground state of Pr $`4f`$-electrons in the crystalline electric field would most likely be the non-Kramers doublet $`\mathrm{\Gamma }_3`$ state. The $`\mathrm{\Gamma }_3`$ ground state has quadrupolar fluctuations and its relations to the heavy fermion behavior and also to a pairing mechanism seems to be possible to discuss. However, its relation to the pairing symmetry, which will be focused in this paper, is unclear at present. Moreover, the singlet $`\mathrm{\Gamma }_1`$ state without quadrupolar fluctuations is also proposed as the ground state. The details of the pairing mechanism will not be discussed here.

The organization of this paper is as follows. In section II, the gap structures for the A- and B-phases are discussed. In section III, jumps of the specific heat at two transition points (Tc1 and Tc2) are calculated. To fit these calculations to the specific heat measurement, one obtains information for the coefficients of the fourth order terms in the GL potential. In section IV, arguments for a spontaneous cubic symmetry breaking in the B-phase are presented, and a summary is given in section V.

## II The Gap structures

In this section, the gap structures for the A- and B-phases in the spin singlet pairing case are discussed within the frameork of the phenomenological GL theory. Information on the gap structure along the $`c`$-axis, which is still unresolved by the thermal conductivity measurement, are obtained.

The GL theory to be considered is based on the cubic crystal symmetry group Th. This group has essentially three irreducible representations A1 ($`\mathrm{\Gamma }_1`$), E ($`\mathrm{\Gamma }_2,\mathrm{\Gamma }_3`$) and T ($`\mathrm{\Gamma }_5`$).

### A A-phase

In this phase, four point-node-like (four-dip) structures of the gap function located within the $`ab`$-plane are observed along the  and  directions . Then, the question is whether the gap function has two additional dips along $`c`$-axis (i.e. the gap has a total of six dips) or not (i.e. four dips in all).

It should be pointed out that Hc2 is isotropic in this superconductor. In general, the superconducting states represented by a multicomponent order parameter should have an anisotropic Hc2 . Therefore, the gap function in the A-phase appears to be represented by a single component order parameter, i.e., the A1 representation.

The gap function with six dips along the ,  and  directions belongs to the A1 representation because this structure has the cubic symmetry Th. Actually, one can see this as follows. Let $`\eta `$ and $`\varphi _s(𝐤)`$ stand for the order parameter and the basis function of the A1 representation, respectively. By the definition, $`\varphi _s(𝐤)`$ should be invariant under Th. The gap function in this representation is

$$\mathrm{\Delta }(𝐤)=\eta \varphi _s(𝐤).$$

The six-dip state can be represented by using invariants:

$$\varphi _s(𝐤)=s+k_x^4+k_y^4+k_z^4(s1),$$
(1)
and also

$$\varphi _s(𝐤)=s+k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2(s<<1),$$
(2)
$`\mathrm{}`$ etc. are possible. When one takes $`s=1`$ in Eq. (1) ($`s=0`$ in Eq. (2)), the dips become point nodes. Gap functions such as in Eqs. (1) and (2) have been discussed in the context of the YNiB2C2 superconductor. It should be noted that the function $`\varphi _s(𝐤)`$ which can give rise to six dips is not limited to these two forms. There are many possible forms which have six dips. The point is that the six-dip state can be described in the A1 representation and essentially categorized into the so-called ”anisotropic $`s`$-wave” state.

The four-dip state does not belong to the A1 representation because this gap structure does not have the cubic symmetry. Actually, one would see in the next subsection that this state belongs to the A$`{}_{1}{}^{}`$E combined representation. Then, a multicomponent orderparameter is needed to describe this gap function (See, Eq. (4)), and in general, an anisotropy in $`H_{c2}`$ is to be expected. in contradiction to the experimental result. Moreover, in the next subsection one would see that the four-dip structure arises in a somewhat accidental way (See, Eq. (5)).

Then, in the A-phase, it seems to be natural to consider that the six-dip state (strongly anisotropic $`s`$-wave state) is realized.

### B B-phase

In this phase, two dips of the gap function within in $`ab`$-plane along  direction are observed. Then, the question is whether the gap function has two additional dips along the $`c`$-axis (i.e. the gap has in total four dips) or not (i.e. a total of two dips).

Phase transitions are observed by varying $`H`$ and also varying $`T`$ without a magnetic field turned on (Fig. I). To describe the phase transition without external fields, one can consider the combined representations. Possible combinations in the present case are A$`{}_{1}{}^{}`$E, A$`{}_{1}{}^{}`$ T, E$``$T, and A$`{}_{1}{}^{}`$ E $``$ T.

As we will see below, both the four-dip and two-dip states belong to the A$`{}_{1}{}^{}`$ E representation. Let $`g_{\mathrm{A}_1}`$ and $`g_\mathrm{E}`$ stand for the effective coupling constant of the pairing interaction for the A1 channel and the E channel, respectively. In general, these values are split by the crystal field. The A$`{}_{1}{}^{}`$ E combined state is realized when

$$g_{\mathrm{A}_1}g_\mathrm{E},$$
(3)
is satisfied. This comes from some accidental degeneracy.

The gap function in A$`{}_{1}{}^{}`$ E model is expanded by the basis functions of A1 and E representations:

$$\mathrm{\Delta }(𝐤)=\eta \varphi _s(𝐤)+\eta _1\varphi _{2y^2z^2x^2}(𝐤)+\eta _2\varphi _{z^2x^2}(𝐤).$$
(4)
The function $`\varphi _{2y^2z^2x^2}(𝐤)`$ is the first component of the 2D basis of the E representation which transforms like the function $`2k_y^2k_z^2k_x^2`$ under the Th group. It has two line nodes parallel to the $`ac`$-plane. The function $`\varphi _{z^2x^2}(𝐤)`$ is the second component of the 2D basis of the E representation which transforms like the function $`k_z^2k_x^2`$ under the Th group. It has two line nodes crossing at two points along the  direction. Note that these functions may include terms beyond quadratic order in the momentum $`𝐤`$. However the nodal structure of these basis functions remain unaltered, in general, even if the higher order terms are included.

To obtain a four-dip state along the  and  directions, one must fine-tune the three basis functions into the particular forms $`\varphi _s(𝐤)=12(k_x^4+k_y^4+k_z^4)/3`$, $`\varphi _{2y^2z^2x^2}(𝐤)=2k_y^4k_z^4k_x^4`$, $`\varphi _{z^2x^2}(𝐤)=\sqrt{3}(k_z^4k_x^4)`$, where the second order terms of $`𝐤`$ are absent, and also the three order parameters are fixed to take the particular values

$$(\eta ,\eta _1,\eta _2)(1,\frac{1}{6},\frac{1}{2\sqrt{3}}).$$
(5)
This looks highly accidental. The gap function with four dips along  and  directions discussed in the previous subsection can also be obtained in the same manner.

Two-dip state is realized when

$$(\eta ,\eta _1,\eta _2)=(\alpha ,0,\pm i\beta ),\alpha ,\beta 𝐑,$$
(6)
because $`\varphi _{z^2x^2}(𝐤)`$ is zero and $`\varphi _s(𝐤)`$ becomes minimum along  direction. Note that the two-dip structure appears generically for all values of the two parameters $`\alpha `$ and $`\beta `$. Moreover, it is not necessary to fix the basal functions $`\varphi _s(𝐤)`$, $`\varphi _{2y^2z^2x^2}(𝐤)`$ and $`\varphi _{z^2x^2}(𝐤)`$ to particular forms. The only assumption needed here is the six-dip structure of the function $`\varphi _s(𝐤)`$. It should be noted again that the function $`\varphi _s(𝐤)`$ is not limited to a particular form to obtain six dips. There are many possible forms to obtain six dips and described in the A1 representation (anisotropic $`s`$-wave state). Eq. (1) (or Eq. (2)) is one of the examples. Then, the state Eq. (6) is an anisotropic $`s+id_{z^2x^2}`$-wave state and breaks the cubic crystal symmetry and also the time reversal symmetry spontaneously.

Let us discuss the stability of these two states Eqs. $`(\text{5})`$ and $`(\text{6})`$. It will be shown that the state Eq. $`(\text{6})`$ is found as a stable state of the system, but the state Eq. $`(\text{5})`$ would be not.

Consider the GL potential for A$`{}_{1}{}^{}`$ E representation of Th

$`F_{\mathrm{op}}`$ $`=`$ $`A_{\mathrm{A}_1}(T)|\eta |^2+\beta |\eta |^4+A_\mathrm{E}(T)(|\eta _1|^2+|\eta _2|^2)`$ (10)
$`+\beta _1(|\eta _1|^2+|\eta _2|^2)^2+\beta _2(\eta _1\eta _2^{}\eta _1^{}\eta _2)^2`$
$`+\theta _1|\eta |^2(|\eta _1|^2+|\eta _2|^2)`$
$`+{\displaystyle \frac{\theta _2}{2}}\{\eta ^2(\eta _1^2+\eta _2^2)+c.c\},`$
where

$$A_{\mathrm{A}_1}(T)=\frac{T}{T_{\mathrm{A}_1}}1,A_\mathrm{E}(T)=\frac{T}{T_\mathrm{E}}1,$$

and

$$T_{\mathrm{A}_1}\mathrm{exp}[1/g_{\mathrm{A}_1}N(0)],T_\mathrm{E}\mathrm{exp}[1/g_\mathrm{E}N(0)],$$

in the weak coupling approximation (See. Eq. (3)). $`N(0)`$ is the density of states at the Fermi surface. We consider a temperature region $`T<min[T_{\mathrm{A}_1},T_\mathrm{E}]`$. The terms proportional to $`\theta _1`$ and $`\theta _2`$ denote the coupling between the order paremeter in the A1 representation and the two order parameters in the E representation. In general, another coupling term

$$\theta _3[e^{i\gamma _3}\eta ^{}(\eta _1|\eta _1|^22\eta _1|\eta _2|^2\eta _1^{}\eta _2^2)+c.c.]$$

could be included in the potential but this term makes the phase transition 1st order, which appears to be inconsistent with the specific heat experiments. So, this term is excluded here. For simplicity, we take

$$\theta _1=\theta _2.$$
(11)
By using the parameterization

$$\eta =|\eta |,\eta _1=|\stackrel{~}{\eta }|\mathrm{cos}\psi e^{i\varphi _1},\eta _2=|\stackrel{~}{\eta }|\mathrm{sin}\psi e^{i\varphi _2},$$

Eq. (10) becomes

$`F_{op}`$ $`=`$ $`A_{\mathrm{A}_1}(T)|\eta |^2+A_\mathrm{E}(T)|\stackrel{~}{\eta }|^2+\beta |\eta |^4+\beta _1|\stackrel{~}{\eta }|^4`$ (14)
$`\beta _2|\stackrel{~}{\eta }|^4\mathrm{sin}^22\psi \mathrm{sin}^2(\varphi _1\varphi _2)+\theta _1|\eta |^2|\stackrel{~}{\eta }|^2`$
$`+\theta _2|\eta |^2|\stackrel{~}{\eta }|^2\left[\mathrm{cos}^2\psi \mathrm{cos}2\varphi _1+\mathrm{sin}^2\psi \mathrm{cos}2\varphi _2\right].`$
After a bit tedious calculation one obtains all of the stable solutions:

1. $`0<\beta _2`$

$`(\eta ,\eta _1,\eta _2)`$ $`=`$ $`(a,{\displaystyle \frac{b}{\sqrt{2}}}e^{i\varphi },\pm i{\displaystyle \frac{b}{\sqrt{2}}}e^{i\varphi }),`$ (15)
where,

$$a=\sqrt{\frac{A_{\mathrm{A}_1}(T)}{2\beta }},b=\pm \sqrt{\frac{A_\mathrm{E}(T)}{2(\beta _1|\beta _2|)}},$$

and $`\varphi `$ is arbitral here.
2. $`\beta _2<0`$

$$(\eta ,\eta _1,\eta _2)=\{\begin{array}{c}(c,\pm id\mathrm{cos}\psi ,\pm id\mathrm{sin}\psi ),\mathrm{for}0<\theta _2\hfill \\ (c,\pm d\mathrm{cos}\psi ,\pm d\mathrm{sin}\psi ),\mathrm{for}\theta _2<0,\hfill \end{array}$$
(16)
where,

$$c=\sqrt{\frac{A_{\mathrm{A}_1}(T)}{2\beta }},d=\sqrt{\frac{A_\mathrm{E}(T)}{2\beta _1}},$$

and $`\psi `$ is arbitral. The degeneracy of $`\psi `$ is lifted by taking into account the 6-th order terms for $`(\eta _1,\eta _2)`$ as a weak perturbation

$`\gamma _1(|\eta _1|^2+|\eta _2|^2)^3+\gamma _2(|\eta _1|^2+|\eta _2|^2)|\eta _1^2+\eta _2^2|^3`$ (17)
$`+\gamma _3|\eta _1|^2|3\eta _2^2\eta _1^2|.`$ (18)
The coupling terms between $`\eta `$ and $`(\eta _1,\eta _2)`$ are neglected here. Then,

$$\psi =\{\begin{array}{c}0,\pi (0<\gamma _3)\hfill \\ \pm \frac{\pi }{2}(\gamma _3<0).\hfill \end{array}$$
(19)

From Eqs. (4), (6), (16) and (19), one can see that the two-dip state (anisotropic $`s+id_{z^2x^2}`$-wave state) is stable when $`\beta _2,\gamma _3<0<\theta _2`$. On the other hand, it is hard to see how the four-dip state Eq. (5) can be stabilized in the present argument.

The gap structures for the other stable solutions (Eq. (15) and Eq. (16) for $`\theta _2<0`$) are inconsistent with the thermal conductivity measurement. Furthermore, based on a consideration for the condensation energy within the weak coupling approach, it has been pointed out that the $`s+id`$-wave state is energetically favored in A$`{}_{1}{}^{}E`$ combined representation.

Fron these considerations, we conclude that the two-dip state appears to be realized in the B-phase.

### C phase transition and symmetry breaking

Let us summarize the results in this section and try to describe the phase transition in the absence of the magnetic field.

Put the parameters

$$T_{\mathrm{A}_1}=T_{c1}=1.85K,T_\mathrm{E}=T_{c2}=1.75K,$$
(20)
and $`\beta _2,\gamma _3<0<\theta _2`$ in the GL free energy Eq. (10). For simplicity, we take Eq. (11). Then, the stationary values of the order parameters in the A-phase ($`T_{c2}<T<T_{c1}`$) are

$$\eta =\sqrt{\frac{A_{\mathrm{A}_1}(T)}{2\beta }},\eta _1=\eta _2=0,$$
(21)
and in the B phase ($`T<T_{c2}`$)

$$\eta =\sqrt{\frac{A_{\mathrm{A}_1}(T)}{2\beta }},\eta _1=0,\eta _2=\pm i\sqrt{\frac{A_\mathrm{E}(T)}{2\beta _1}}.$$
(22)
From Eqs. (4), (21) and (22), one can see that the gap function in the A-phase has six dips along the ,  and  directions (anisotropic $`s`$-wave state), and that in the B-phase there are two dips along the  direction (anisotropic $`s+id_{z^2x^2}`$-wave state). The result is consistent with the thermal conductivity measurement.

The gap function in the B-phase breaks the cubic crystal symmetry and also the time reversal symmetry. The consequence of this cubic symmetry breaking will be discussed in section IV. As a signal of the time reversal symmetry breaking, it is expected that a spontaneous magnetization be observed by the $`\mu `$SR measurement .

## III Specific heat jump

In this section, we calculate the jumps of the specific heat at the two transition points $`T_{c1}=1.85K`$ and $`T_{c2}=1.75K`$ in Fig. I indicated by the specific heat measurement. To fit the calculation to the experimental result the ratio of the coefficients of 4-th order terms in GL potential $`\beta /\beta _1`$ is obtained. The simplification Eq. (11) is used.

Let us introduce

$$\stackrel{}{\eta }=(\eta ,\eta _1,\eta _2).$$

By using the stationary condition for the free energy Eq. (10)

$$0=\frac{F_{op}}{\stackrel{}{\eta }},$$
(23)
the specific heat is

$`C(T)`$ $`=`$ $`C_0(T)T{\displaystyle \frac{^2F_{op}}{T^2}}`$ (24)
$`=`$ $`C_0(T)T\{{\displaystyle \frac{\stackrel{}{\eta }}{T}}{\displaystyle \frac{^2F}{\stackrel{}{\eta }T}}+c.c.\},`$ (25)
where $`C_0(T)`$ is the background specific heat. Let us define $`\mathrm{\Delta }C(T)=C(T)C_0(T)`$. Then, in the A-phase,

$$\mathrm{\Delta }_\mathrm{A}C(T)=\frac{T}{2\beta T_1^2},$$
(26)
while in the B-phase,

$$\mathrm{\Delta }_\mathrm{B}C(T)=T\left(\frac{1}{2\beta T_1^2}+\frac{1}{2\beta _1T_3^2}\right).$$
(27)
Therefore, the ratio of the jumps at the two transition points is

$$\frac{\mathrm{\Delta }_\mathrm{B}C(T_{c2})\mathrm{\Delta }_\mathrm{A}C(T_{c2})}{\mathrm{\Delta }_\mathrm{A}C(T_{c1})}=\frac{\beta T_{c1}}{\beta _1T_{c2}}$$
(28)
The experimental result for this quantity is around 1 . Then, from Eq. (20), one can obtain the ratio

$$\beta /\beta _10.9.$$
(29)
## IV Cubic symmetry breaking in B-phase

the B-phase gap breaks the cubic crystal symmetry $`T_h`$ spontaneously. So, one can expect a crystal deformation with the cubic symmetry breaking. In this section, the strain-order parameter coupling is taken into account and thermal expansions and isothermal elastic constants are calculated. In the B-phase, one obtains anisotropic results for these values. This comes from the spontaneous crystal symmetry breaking.

### A strain-orderparameter coupling

The GL potential which include the strain-order parameter coupling is,

$`F_{\mathrm{GL}}`$ $`=`$ $`F_{\mathrm{op}}+F_{\mathrm{st}\mathrm{op}}+F_{\mathrm{el}},`$ (30)
$`F_{\mathrm{st}\mathrm{op}}`$ $`=`$ $`|\eta |^2C_{\mathrm{A}_1}(ϵ_{xx}+ϵ_{yy}+ϵ_{zz})`$ (38)
$`|\eta _1|^2[C_{\mathrm{A}_1}(ϵ_{xx}+ϵ_{yy}+ϵ_{zz})`$
$`+C_\mathrm{E}(2ϵ_{yy}ϵ_{zz}ϵ_{xx})]`$
$`|\eta _2|^2[C_{\mathrm{A}_1}(ϵ_{xx}+ϵ_{yy}+ϵ_{zz})`$
$`C_\mathrm{E}(2ϵ_{yy}ϵ_{zz}ϵ_{xx})]`$
$`\sqrt{3}C_\mathrm{E}(\eta _1\eta _2^{}+c.c.)(ϵ_{zz}ϵ_{xx})`$
$`C_{\mathrm{AE}}(\eta \eta _1^{}+c.c.)(2ϵ_{yy}ϵ_{zz}ϵ_{xx})`$
$`\sqrt{3}C_{\mathrm{AE}}(\eta \eta _2^{}+c.c.)(ϵ_{zz}ϵ_{xx}),`$
$`F_{\mathrm{el}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ijkl}{}}c_{ijkl}^{(0)}ϵ_{ij}ϵ_{kl},`$ (39)
where $`c_{ijkl}^{(0)}`$ is the background elastic constants. Obviously, it is symmetric under the Th symmetry.

Let us consider the amplitude of the strain. It is obtained by solving the equation

$$0=\frac{F_{GL}}{ϵ_{ij}},$$

We introduce the notation $`ϵ_{xx}=ϵ_1,ϵ_{yy}=ϵ_2,ϵ_{zz}=ϵ_3,ϵ_{xy}=ϵ_4,ϵ_{yz}=ϵ_5,ϵ_{zx}=ϵ_6.`$ In the A-phase (see, Eq. (21)), one finds the isotropic result

$`ϵ_1`$ $`=`$ $`ϵ_2=ϵ_3={\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{A_{\mathrm{A}_1}(T)}{2\beta }},`$ (40)
$`ϵ_4`$ $`=`$ $`ϵ_5=ϵ_6=0,`$ (41)
where $`\lambda =c_{11}^{(0)}=c_{22}^{(0)}=c_{33}^{(0)}`$. In B-phase (See, Eq. (22)), an anisotropy can be seen in the $`y`$ direction:

$`ϵ_3`$ $`=`$ $`ϵ_1={\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{A_{\mathrm{A}_1}(T)}{2\beta }}{\displaystyle \frac{C_{\mathrm{A}_1}+C_E}{\lambda }}{\displaystyle \frac{A_\mathrm{E}(T)}{2\beta _1}},`$ (42)
$`ϵ_2`$ $`=`$ $`{\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{A_{\mathrm{A}_1}(T)}{2\beta }}{\displaystyle \frac{C_{\mathrm{A}_1}2C_E}{\lambda }}{\displaystyle \frac{A_\mathrm{E}(T)}{2\beta _1}},`$ (43)
$`ϵ_4`$ $`=`$ $`ϵ_5=ϵ_6=0.`$ (44)
It should be noted that the terms proportional to $`C_{\mathrm{AE}}`$ in Eq. (39) play no roles in the state Eq. (22).

### B thermal expansion

The presence of the strain can be measured by the thermal expansion which is definition by

$$\alpha _\mu =\alpha _\mu ^{(0)}+\frac{ϵ_\mu }{T},$$
(45)
where $`\mu =1,\mathrm{},6`$, and $`\alpha _\mu ^{(0)}`$ is the background value. Define $`\mathrm{\Delta }\alpha _\mu =\alpha _\mu \alpha _\mu ^{(0)}`$. Then, in the A-phase (See, Eq. (41)),

$`\mathrm{\Delta }_\mathrm{A}\alpha _1`$ $`=`$ $`\mathrm{\Delta }_\mathrm{A}\alpha _2=\mathrm{\Delta }_\mathrm{A}\alpha _3={\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{1}{2\beta T_{c1}}}`$ (46)
$`\mathrm{\Delta }_\mathrm{A}\alpha _4`$ $`=`$ $`\mathrm{\Delta }_\mathrm{A}\alpha _5=\mathrm{\Delta }_\mathrm{A}\alpha _6=0.`$ (47)
In B-phase (See, Eq. (44)), the thermal expansion becomes anisotropic:

$`\mathrm{\Delta }_\mathrm{B}\alpha _3`$ $`=`$ $`\mathrm{\Delta }_\mathrm{B}\alpha _1={\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{1}{2\beta T_{c1}}}{\displaystyle \frac{C_{\mathrm{A}_1}+C_\mathrm{E}}{\lambda }}{\displaystyle \frac{1}{2\beta _1T_{c2}}},`$ (48)
$`\mathrm{\Delta }_\mathrm{B}\alpha _2`$ $`=`$ $`{\displaystyle \frac{C_{\mathrm{A}_1}}{\lambda }}{\displaystyle \frac{1}{2\beta T_{c1}}}{\displaystyle \frac{C_{\mathrm{A}_1}2C_\mathrm{E}}{\lambda }}{\displaystyle \frac{1}{2\beta _1T_{c2}}},`$ (49)
$`\mathrm{\Delta }_\mathrm{B}\alpha _4`$ $`=`$ $`\mathrm{\Delta }_\mathrm{B}\alpha _5=\mathrm{\Delta }_\mathrm{B}\alpha _6=0.`$ (50)
The thermal expansion along  direction has been measured, and two clear jumps are observed at $`T_{c1}`$ and $`T_{c2}`$. The ratio of these two jumps is

$$\frac{\mathrm{\Delta }_\mathrm{B}\alpha _1\mathrm{\Delta }_\mathrm{A}\alpha _1}{\mathrm{\Delta }_\mathrm{A}\alpha _1}=\left(1+\frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}\right)\frac{\beta }{\beta _1}\frac{T_{c1}}{T_{c2}}1.8.$$
(51)
From Eqs. (20) and (29), one can obtain a relation

$$C_\mathrm{E}/C_{\mathrm{A}_1}0.9,$$
(52)
and estimate the ratio of jumps at the two transition points for the thermal expansion in the  direction

$$\frac{\mathrm{\Delta }_\mathrm{B}\alpha _2\mathrm{\Delta }_\mathrm{A}\alpha _2}{\mathrm{\Delta }_\mathrm{A}\alpha _2}=\left(12\frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}\right)\frac{\beta }{\beta _1}\frac{T_{c1}}{T_{c2}}0.8.$$
(53)
This quantity has not been measured yet.

### C isothermal elastic constants

The isothermal elastic constant is

$$c_{\mu \nu }=\left(\frac{^2F_{GL}}{ϵ_\mu ϵ_\nu }\right)_T.$$
(54)
The small variation of the strain changes the amplitude of the order parameters as

$`\delta |\eta |^2`$ $`=`$ $`{\displaystyle \frac{C_{\mathrm{A}_1}}{2\beta }}(\delta ϵ_1+\delta ϵ_2+\delta ϵ_3),`$ (55)
$`\delta |\eta _2|^2`$ $`=`$ $`{\displaystyle \frac{C_{\mathrm{A}_1}}{2\beta _1}}(\delta ϵ_1+\delta ϵ_2+\delta ϵ_3)`$ (57)
$`{\displaystyle \frac{C_\mathrm{E}}{2\beta _1}}(2\delta ϵ_2\delta ϵ_3\delta ϵ_1).`$
By using the condition Eq. (23),

$$c_{\mu \nu }=c_{\mu \nu }^{(0)}+[\frac{\stackrel{}{\eta }}{ϵ_\mu }\frac{^2F_{\mathrm{GL}}}{\stackrel{}{\eta }ϵ_\nu }+c.c.],$$
(58)
where $`c_{\mu \nu }^{(0)}`$ is the background elastic constant. Let us define $`\mathrm{\Delta }c_{\mu \nu }=c_{\mu \nu }c_{\mu \nu }^{(0)}`$. In A-phase,

$`\mathrm{\Delta }_\mathrm{A}c_{11}=\mathrm{\Delta }_\mathrm{A}c_{22}=\mathrm{}=\mathrm{\Delta }_\mathrm{A}c_{31}={\displaystyle \frac{C_\mathrm{A}^2}{2\beta }},`$ (59)
and the others are zero. In B-phase,

$`\mathrm{\Delta }_\mathrm{B}c_{33}`$ $`=`$ $`\mathrm{\Delta }_\mathrm{B}c_{11}=\mathrm{\Delta }_\mathrm{B}c_{31}={\displaystyle \frac{C_{\mathrm{A}_1}^2}{2\beta }}{\displaystyle \frac{(C_{\mathrm{A}_1}+C_\mathrm{E})^2}{2\beta _1}},`$ (60)
$`\mathrm{\Delta }_\mathrm{B}c_{22}`$ $`=`$ $`{\displaystyle \frac{C_{\mathrm{A}_1}^2}{2\beta }}{\displaystyle \frac{(C_{\mathrm{A}_1}2C_\mathrm{E})^2}{2\beta _1}},`$ (61)
$`\mathrm{\Delta }_\mathrm{B}c_{23}`$ $`=`$ $`\mathrm{\Delta }_\mathrm{B}c_{12}={\displaystyle \frac{C_{\mathrm{A}_1}^2}{2\beta }}{\displaystyle \frac{(C_{\mathrm{A}_1}+C_\mathrm{E})(C_{\mathrm{A}_1}2C_\mathrm{E})}{2\beta _1}},`$ (62)
and the others are zero. From Eqs. (29) and (52), one can estimate the ratios of jumps at two transition points of the elastic constants,

$`{\displaystyle \frac{\mathrm{\Delta }_\mathrm{B}c_{11}\mathrm{\Delta }_\mathrm{A}c_{11}}{\mathrm{\Delta }_\mathrm{A}c_{11}}}`$ $`=`$ $`{\displaystyle \frac{\beta }{\beta _1}}\left(1+{\displaystyle \frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}}\right)^23.6,`$ (63)
$`{\displaystyle \frac{\mathrm{\Delta }_\mathrm{B}c_{22}\mathrm{\Delta }_\mathrm{A}c_{22}}{\mathrm{\Delta }_\mathrm{A}c_{22}}}`$ $`=`$ $`{\displaystyle \frac{\beta }{\beta _1}}\left(12{\displaystyle \frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}}\right)^20.6,`$ (64)
$`{\displaystyle \frac{\mathrm{\Delta }_\mathrm{B}c_{12}\mathrm{\Delta }_\mathrm{A}c_{12}}{\mathrm{\Delta }_\mathrm{A}c_{12}}}`$ $`=`$ $`{\displaystyle \frac{\beta }{\beta _1}}\left(1+{\displaystyle \frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}}\right)\left(12{\displaystyle \frac{C_\mathrm{E}}{C_{\mathrm{A}_1}}}\right).`$ (65)
$``$ $`1.4.`$ (66)
In conclusion for this section, the thermal expansion and the elastic constant are isotropic in the A-phase, but anisotropic in the B-phase. To observe the anisotropy would serve as additional evidences for the phase transition and the spontaneous cubic symmetry breaking in the B-phase.

## V Summary and discussions

In summary, the multiple superconducting phases (A-phase and B-phase) in the skutterudite PrOs4Sb12 with cubic crystal symmetry Th is discussed by using the phenomenological Ginzburg-Landau approach. The spin singlet pairing case is considered here. According to the thermal conductivity measurement , the gap function has within in the $`ab`$-plane four point-node-like structure (four dips) along the  and  directions on the Fermi surface in the A-phase, and two dips along the  direction in the B-phase. But the gap structure along the $`c`$-axis was unresolved. This problem is considered here.

Because of the absence of an anisotropy in Hc2, the superconductivity in A-phase appears to be the conventional state (anisotropic $`s`$-wave state). Then, the gap structure has the cubic symmetry (the A1 representation) and leads to two additional dips along the $`c`$-axis. In the B-phase, a state without additional dips along the $`c`$-axis (anisotropic $`s+id_{z^2x^2}`$-wave state) is stabilized by the Ginzburg-Landau potential of the A1$``$E combined representation in a wide parameter region. On the other hand, it is hard to see that a state with two dips along $`c`$-axis is stabilized, since to obtain such a state requires a highly accidental mixing between the basis functions of the crystal symmetry Th. Moreover, it has been shown that the $`s+id`$-wave state is energetically favored in the A$`{}_{1}{}^{}`$E combined representation based on the consideration for the condensation energy of each state within the weak coupling approach. Hence it is natural to expect that the gap function in A-phase has a total of six dips along , and  directions and while in the B-phase there are in total two dips along the  direction. The phase transition can be described by the Ginzburg-Landau potential.

The gap function inthe B-phase breaks the cubic crystal symmetry Th spontaneously. This symmetry breaking would cause anisotropic crystal strain and anisotropic anomaly of the thermal expansion and the isothermal elastic constant could be observed .

The gap function in B-phase breaks the time reversal symmetry (See, Eqs. (4) and (6)). Then, the spontaneous magnetization around impurities is expected to be observed by the $`\mu `$SR measurement. The magnetization have not been detected yet. More accurate measurements are highly desired.

Another expected phenomenon in the B-phase is anomalous flux flow which comes from the domain formation, because the domain structure would exist in the B-phase since the stable state Eq. (6) is two-fold degenerate.

Additional phase transitions induced by uniaxial pressure may be possible in this superconductor. Such a transition could occur in a superconductor described by multicomponent order parameter.

###### Acknowledgements.

The author is grateful to N. Hayashi, K. Izawa, K. Maki, Y. Matsuda, Y. Nakajima, N. Oeschler, M. Sigrist, A. Tanaka, K. Ueda, I. Vekhter, Q. Yuan, and especially to P. Thalmeier for useful discussions and comments.


text:
# A cellular automaton model of gravitational clustering

\[

## Abstract

Gravitational clustering of a random distribution of point masses is dominated by the effective short-range interactions due to large-scale isotropy. We introduce a one-dimensional cellular automaton to reproduce this effect in the most schematic way: at each time particles move towards their nearest neighbours with whom they coalesce on collision. This model shows an extremely rich phenomenology with features of scale-invariant dynamics leading to a tree-like structure in space-time whose topological self-similarity are characterised with universal exponents. Our model suggests a simple interpretation of the non-analytic hierarchical clustering and can reproduce some of the self-similar features of gravitational N-body simulations.

PACS numbers:89.75.Hc, 61.43.Hv, 98.80.-k

\]

The basic mechanism of large-scale pattern formation by gravity remains mainly elusive, in spite of elaborate studies in the past few decades . Models of pattern formation, such as diffusion-limited aggregation (DLA) and its many variants , which have been remarkably successful in wide ranges of disciplines, do not appear to represent, even at a simplified level, gravitational clustering. A different class of irreversible aggregation models, described by Smoluchowski coagulation equation and its various associated kernels, gives rise to a power-law distribution of the masses of the aggregates which are, however, randomly distributed in space .

Gravitational clustering, in general, leads to a power law distribution of the masses, known as the Press-Schechter mass function , but with a space distribution which, at least up to some scale, is not homogeneous and can be described by fractal geometry . The key difference between the aggregations one encounters in statistical physics and gravitational clustering is that in the former, forces are short-ranged, and the nonlinear dynamics, driven by collisions, erase the memory of the initial conditions rather fast. Gravity, on the other hand, is long-ranged with a deterministic evolution which strongly depends on the initial conditions even after the dynamics become nonlinear.

Despite these notorious features of gravity, it has long been shown by Chandrasekhar  that for the Poisson distribution of particle positions, gravity is effectively short-ranged : it can be well-described by a simple nearest-neighbours interactions. For the Poisson distribution of particles, the magnitude of the gravitational force, $`|𝐅|`$, has a Holtsmark distribution  given by

$`W(|𝐅|)`$ $`=`$ $`{\displaystyle \frac{2|𝐅|}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑yy\mathrm{sin}\left(y|𝐅|\right)e^{ay^{3/2}}`$ (1)
$``$ $`{\displaystyle \frac{1}{|𝐅|^{5/2}}}\mathrm{as}|𝐅|\mathrm{}`$ (2)
where $`a=4n(2\pi Gm)^{3/2}/15`$, which is a constant depending on the mass of each particle, $`m`$, and the uniform number density $`n`$ . It was shown analytically and verified numerically that the $`5/2`$ power-law tail of the Holtsmark distribution, is exactly due to the forces between the nearest neighbours . Even more remarkable is the fact that the distribution of the forces from the first neighbour approximation agrees with (2) over most of the range of $`|𝐅|`$ . That the long-range gravitational force can be well-described by an effective short-range interaction may at first appear puzzling. However, the reason for this simplification is that the isotropy of Poisson distribution is only broken at small scales where the granularity becomes important.

As the system evolves, particle positions become correlated, but the self-similar nature of gravitational clustering (see for example ) can ensure that at progressively larger scales the functional form of the Holtsmark distribution remains intact. Hence, under an appropriate renormalisation of the masses and the distances, determined by the scale of granularity at a given time, gravitational interactions can be taken to be effectively short-ranged, even at times comparable to the dynamical time of the system.

Inspired by these facts, we introduce an extremely simple one-dimensional automaton model which focuses precisely on this granular and self-similar nature of the gravitational clustering phenomenon. Unlike most previous models such as DLA, Smoluchowski or 1-dimensional Burgers, our aggregation rule is position-dependent and independent of the masses and initial velocities of the aggregates. Unlike the conventional boolean cellular automata, we distribute particles randomly on a line rather than on discrete lattice points. We focus mainly on the universal features of the distribution of the masses of the aggregates and on the space-time topological properties of the full aggregation history, i.e., on the universal characteristics of the “merger tree” of our model.

The basic idea of this model, that the particles move towards their nearest neighbours, is implemented by the following algorithm. We distribute $`10^5`$ particles randomly on a periodic line with a length of $`10^5`$ units, so that the average density is equal to one. Particles move towards their nearest neighbours either by one unit at each time step or by half their separations, whichever that is shorter. If two particles are closer than a lower threshold, which we set equal to 2 units, and, in addition, are mutual nearest-neighbours they coalesce at the mid-point of their separation, conserving mass. This aggregation rule for a particle $`i`$ at the position $`x_i(t)`$ at time $`t`$ can be written as,

$`x_i(t+\mathrm{\Delta }t)`$ $`=`$ $`x_i(t)+\mathrm{Min}[{\displaystyle \frac{x_r(t)}{2}},1]\mathrm{If}x_r(t)<|x_l(t)|`$ (3)
$`x_i(t+\mathrm{\Delta }t)`$ $`=`$ $`x_i(t)+\mathrm{Max}[{\displaystyle \frac{x_l(t)}{2}},1]\mathrm{If}x_r(t)>|x_l(t)|`$ (4)
where $`x_r(t)=x_{i+1}(t)x_i(t)`$ and $`x_l(t)=x_{i1}(t)x_i(t)`$. This process is illustrated in Fig. 1.

As time progresses, masses and distances of the nearest neighbours rescale and the same rule, given by (LABEL:rule), is followed by more massive and farther apart aggregates. The continuation of this process will eventually lead to the total coagulation of the colloidal particles into one single mass. This aggregation mechanism traces a self-similar tree-like structure in space-time as shown in Fig. 2.

The tree structure of the aggregation process in space-time is a manifestation of topological self-similarity , which is a property of many branched structures such as river-networks  and bronchial trees  and can be quantified by various scaling exponents. A most common of these exponents is the Strahler index, $`s`$, given by

$$P(\eta )\eta ^s,$$
(6)
where $`P(\eta )`$ is the probability that a point in the network is connected to $`\eta `$ other points uphill, also known as the drainage area. The Strahler index, $`s`$, is a measure of the bushiness of a branched structure and has an upper limit of $`1.5`$ for river networks .

In our model, the function $`P(\eta )`$, given by (6), interprets as the probability that a newly-formed aggregate has a mass equal to $`\eta `$. We observe a larger value for Strahler index ($`s=2`$, given by the slope of the plot in the upper inset of Fig. 3) than is expected for river networks and which seems to be a characteristic of Cayley tree structures. A notable example of such a structure, also with exponent $`2`$, has recently been observed for the property of internet connections. The difference between our model and the conventional river-networks could be the stringent requirement of mass conservation here.

Topological scaling, is believed to emerge from a self-similar growth process. Self-similar growth, or dynamical scaling, is a dominant feature of hierarchical gravitational clustering , and is the fundamental reason for the emergence of distribution functions such as Press-Schechter mass function in cosmology. An appropriate way to analyse dynamical scaling is, indeed, to study the mass distribution function $`n(m,t)`$. Note that our function differs from the usual number density by a constant factor which is given by the size of our system.

It can be easily inferred from Fig. 2, and we have verified numerically that the average mass, $`m`$, and the mass variance, $`(m^2m^2)^{1/2}`$ grow linearly with time. It is worth comparing this with the growth of average mass in one-dimensional Burgers, for example with uniform initial velocities and positions, where an exponent of $`2/3`$ has been obtained . The linear growth of average mass also holds for Smoluchowski equation with a constant kernel . The universality between our model and Smoluchowski arises inspite of the fact that particle trajectories are not Brownian here. In addition to this common feature, our model and Smoluchowski equation have similar asymptotic states. In both of these models, clusters collide until all the mass falls into one final clump, whereas in 1-dimensional Burgers, the asymptotic state can contain many clumps with zero momentum.

The distribution of the masses deviates slightly from a simple Gaussian for small masses where it develops a power-law as is shown in Fig. 4. The inset clearly demonstrates that a self-similar condensation process, rather similar to the one observed in gravitational N-body simulations (compare with Fig. 1 of ), has set in. In the process of coagulation, the shape of the mass spectrum seems to become fixed, and the curves move in parallel to the right (increasing aggregate mass).

The results presented so far provide strong evidence that the mass distribution function has the general scaling form:

$$n(m,t)m^\alpha t^\beta \mathrm{exp}\left(m^\gamma /t^\lambda \right),$$
(7)
where the value of the exponents $`\alpha ,\beta ,\gamma `$ and $`\lambda `$ will be found in what follows.

The conservation of the total mass in our model leads to the exponent identity

$$\beta +\frac{\lambda }{\gamma }\left(2+\alpha \right)=0.$$
(8)
In addition, the linear rate of growth of the average mass, $`mt`$, leads to the further scaling identity,

$$\frac{\lambda }{\gamma }=1.$$
(9)
At this point one can already see the emergence of a topological scaling, namely a power-law behaviour in the time-integrated mass distribution function, at small masses. We comment that this is different from Strahler index we found in Fig. 3, since the latter refers only to newly-formed aggregates. This topological scaling, namely the power-law behaviour at small masses of the time-integrated mass distribution function, $`n(m,t)𝑑t1/m^l,`$ is implied by the scaling relation (8), which fixes the value of the exponent $`l`$ to $`l=1`$. This value is also confirmed by our numerical simulations.

A further test of our scaling identities (8) and (9), is provided by the evolution of the maxima, $`n_{\mathrm{max}}`$, of the mass distribution function, i.e., by the rate of decay of the peaks of the main plot in Fig. 4. Using the fact that the average mass grows linearly with time, in the scaling expression (7), we obtain $`n_{\mathrm{max}}(m,t)1/t^2,`$ which is again confirmed by our simulations.

Thus, the scaling identities (8) and (9), reduce our scaling ansatz (7) to the self-similar form

$$t^2n(m,t)\left(\frac{m}{t}\right)^\alpha e^{(m/t)^\gamma }.$$
(10)
The factor of $`1/t^2`$ in (10) is a consequence of mass conservation and the linear growth rate of the average mass, and has also been observed for the constant-kernel solution to Smoluchowski equation : $`n(m,t)=4t^2\mathrm{exp}(2m/t)`$. This solution of Smoluchowski equation, is obtained from our general solution (10) by the transformation $`t2t`$ and by using the following values of the exponents: $`\alpha =0`$, $`\gamma =1`$. We shall soon show that these exponents take different values in our model.

To put our notation in accordance with that used in cosmology, we replace the time variable by the cut-off mass $`m^{}(t)`$. As we have mentioned previously, our simulations show that a typical mass grows linearly with time, i.e. $`m^{}t`$. The exponents $`\alpha `$ and $`\gamma `$ in expression (10) can be found by plotting $`n(m,t)m_{}^{}{}_{}{}^{2}`$ against the ratio $`m/m^{}`$, which we have done in Fig. 5.

We show in Fig. 5 that for over two decades in time scale the functional form of the mass function, given by the RHS of (10), is preserved. The fit in Fig. 5 sets the value of the unknown exponents in (10) to $`\alpha 3`$ and $`\gamma 2`$ and we finally arrive at the approximate expression

$$n(m,t)\frac{1}{m^{}(t)^2}\left(\frac{m}{m^{}(t)}\right)^3e^{(m/m^{}(t))^2}$$
(11)
for the mass distribution function. We emphasis that, our solution is different from the constant-kernel solution to Smoluchowski equation which is an special case of our general scaling solution (10), as previously noted. We comment that unlike some of the previous aggregation models used in cosmology , our distribution function cannot be formed from a white noise initial spectrum. It remains to be seen if, as for Smoluchowski with additive kernel, the introduction of a mass-dependent factor in our aggregation rule (LABEL:rule), would give rise to a Press-Schechter type mass function.

We have also analysed the density-density correlation function which has a power-law behaviour with exponent $`1`$ at small scales, indicating that the mass is distributed on zero-dimensional objects, and a crossover to a constant value at large scales where it reminisces the initial Poisson distribution. The crossover length increases linearly with time and the growth of the structures is dominated by the granular properties which are shifted from small to large scales. Thus, our model differs from the usual statistical models which generate fractals, where large-scale structures are built while substructures are preserved and not destroyed as is the case here. In this sense the present model does not generate asymptotic fractal structures in space.

In conclusion, the seminal result of Chandrasekhar, that the long-range gravitational interactions between randomly distributed particles can be almost exactly replaced by nearest-neighbour interactions, stimulated us to present a simple aggregation model which captures this profound feature of gravitating systems. We have shown that our model exhibits topological self-similarity over many decades of mass scale and dynamical scaling over many decades of temporal scale. These properties make it a simple and intuitive model for the study of gravitational hierarchical clustering.

We thank J. Bertoin, G. Caldarelli, U. Frisch, N. Menci, G. Murante, R. Sheth and A. Stella for useful discussions and suggestions. R. M. was supported by TMR network “Fractal structures and Self-organization” under the contract FMRXCT980183.