text:
# Seismic study of stellar convective cores

## 1 Introduction

It has been shown that a discontinuity in the derivatives of the sound speed in the stellar interior, like those arising at the edge of a convective region, introduces a characteristic oscillatory signature in the frequencies of low degree modes (Gough 1990; Monteiro et al. 1994; Roxburgh & Vorontsov 1994; Basu et al. 1994). From this oscillatory signature it has been possible to put limits on the extent of overshoot below the solar convection zone (Basu 1997). Monteiro et al. (2000) have pointed out that this oscillatory signature can be used to study the location of the base of the convection zone as well as the extent of overshoot below this base in other stars using asteroseismic data for only low degree modes. Monteiro et al. (1998) have also suggested that the same oscillatory signal can be used to study the size of convective cores in massive stars. However, it is not clear if such a signal has been demonstrated using actual stellar models with convective cores.

Similarly, helioseismic inversions for the rotation rate in the solar interior (Thompson et al. 1996; Schou et al. 1998) have shown the presence of a tachocline (Spiegel & Zahn 1992) around the base of the convection zone, where the rotation rate varies from differential rotation inside the convection zone to almost solid body like rotation in the radiative interior. This sharp transition in the rotation rate introduces an oscillatory signature in the frequency splitting coefficients (Mazumdar & Antia 2001), which can be used to study the characteristics of tachoclines in other stars.

In this work we attempt to study the possible oscillatory signal arising from convective cores in stars more massive than the Sun. However, we do not find any significant oscillatory signature of the core and try to identify the reason for the apparent absence of this signal. We also study how the amplitude of oscillatory signal from the outer convection zone varies with stellar mass and age. Since even the signal from the outer convection zone is rather weak, such studies may help us in identifying promising stars where attempts can be made to detect such an oscillatory signal. In the absence of a substantial oscillatory signal from the convective core of stars, we try to identify alternate means to detect the presence of a convective core and to measure its size.

## 2 The technique

In the asymptotic limit (Christensen-Dalsgaard & Berthomieu 1991) the eigenfunctions of stellar oscillation modes have the general form $`\mathrm{cos}(\omega \tau +\varphi )`$, where $`\tau `$ is the acoustic depth given by

$$\tau =_r^R\frac{1}{c(r)}dr.$$
(1)
Here $`R`$ is the stellar radius. Using this form for the eigenfunctions it can be shown (Monteiro et al. 1994) that if we have a discontinuity in the derivatives of the sound speed at a point corresponding to an acoustic depth $`\tau _c`$ then the frequencies of p-modes can be written as

$$\omega _{n,\mathrm{}}=\omega _{n,\mathrm{}}^{(s)}+A(\omega _{n,\mathrm{}})\mathrm{sin}(2\omega _{n,\mathrm{}}\tau _c+\varphi ),$$
(2)
where $`n`$ is the radial order and $`\mathrm{}`$ the degree of the mode and $`\omega _{n,\mathrm{}}^{(s)}`$ is the smooth part of the frequency arising from the smooth variation in the sound speed with depth, while the second term is the oscillatory component. The amplitude $`A(\omega )`$ is a smooth function of $`\omega `$. Here we have neglected terms involving the degree $`\mathrm{}`$ since their effects are small for low degree ($`\mathrm{}3`$) modes. In this work we are primarily interested in stars other than the Sun, where it is not possible to detect oscillations with higher degree. Similar oscillatory terms can arise in the splitting coefficients if there is any steep variation in the rotation rate such as what happens in the tachocline region (Mazumdar & Antia 2001).

Following Basu et al. (1994), we take the fourth difference of the frequencies with respect to $`n`$ to enhance the oscillatory signal. Another advantage of taking the fourth difference is that the smooth part of the frequencies becomes negligible and we do not need to include it in our analysis. This will of course, depend on the smooth component of variation of the sound speed, but at least for the Sun this component is found to be negligible. The sharp variation in sound speed inside the second helium ionisation zone also gives rise to an oscillatory signal, but in this case the region with reduced adiabatic index, $`\mathrm{\Gamma }_1`$, has a small but finite width and in some sense it can be considered as two steps of opposite sign located close to each other. Monteiro & Thompson (1998) have shown that the amplitude of such a signal would be modulated with an oscillatory factor of the form $`\mathrm{sin}^2(\omega \beta )`$, where $`\beta `$ is the acoustic half-width of the He ii ionisation zone. We have not included this factor in our adopted form for the oscillatory signal. Thus for the Sun there are two oscillatory components – one arising from the base of the convection zone and the other from the second helium ionisation zone. A similar form may be expected for other stars with $`M1.2M_{}`$, where there is no convection in the core. For more massive stars with convective cores we may expect a third component arising from the boundary of the convective core.

Since there are no observations yet of stellar oscillations with sufficient accuracy to detect this oscillatory signal, in this work we construct stellar models with different mass and age and calculate their oscillation frequencies. Using the computed frequencies we attempt to check if the location of these regions of rapid variation in sound speed can be determined using the frequencies of only low degree modes, which are likely to be observed by forthcoming asteroseismic missions.

Apart from frequencies, we also study the oscillatory signal in the splitting coefficients arising from possible tachoclines in stars. There is no theory to determine the position of tachoclines in stars, but for simplicity we generally assume that the tachocline is located at the base of the outer convection zone, which is in fact the case for the Sun (Spiegel & Zahn 1992). However, in some cases we have studied the effect of varying the position of the tachocline to see how the oscillatory signal varies with its depth. The main advantage of using the splitting coefficients arising from an assumed tachocline is that, in this case there is only one layer where the rotation rate changes rapidly and its location can be altered to study how the signal varies. This exercise is not possible with frequencies where there are multiple layers with rapid variation in the sound speed and their locations cannot be easily shifted over the entire stellar radius.

In order to study the oscillatory signature in splitting coefficients we assume a model tachocline rotation profile of the form (Mazumdar & Antia 2001)

$$\mathrm{\Omega }(r,\theta )=\frac{\delta \mathrm{\Omega }(5\mathrm{cos}^2\theta 1)}{1+\mathrm{exp}[(r_dr)/w]},$$
(3)
where $`\delta \mathrm{\Omega }`$ is the extent of variation in the rotation rate across the tachocline, lying at a radial position of $`r_d`$ and having a half-width of $`w`$. Using this model rotation profile we can calculate the corresponding splitting coefficients for a star. Because of the choice of latitudinal dependence only the splitting coefficient $`c_3(n,\mathrm{})`$ is found to be non-zero. Neglecting the smooth part, the fourth difference can be written as

$$\delta ^4c_3(n,\mathrm{})=\left(a_0+\frac{a_1}{\omega _{n,\mathrm{}}}+\frac{a_2}{\omega _{n,\mathrm{}}^2}\right)\mathrm{sin}(2\omega _{n,\mathrm{}}\tau +\varphi ).$$
(4)
The parameters $`a_0,a_1,a_2,\tau `$ and $`\varphi `$ can be determined by a nonlinear least squares fit. For the frequencies also we use the same form for fitting except that there are multiple oscillatory terms. Each oscillatory part is defined by a similar set of parameters, but the number of parameters to be fitted is larger and hence it is more difficult to obtain a proper fit. For distant stars it will be possible to detect modes with $`\mathrm{}=0,1,2,3`$ only and hence we use only these modes in our study. As mentioned earlier, the amplitude of the signal arising from the He ii ionisation zone has an oscillatory factor and may not be properly represented by the form we have adopted. However, the half-width $`\beta `$ being quite small, the ‘wavelength’ of oscillation in amplitude would be quite large and in the limited frequency range that is likely to be observed we expect less than half a wavelength to be accommodated. Such a variation can be approximated reasonably well by the three terms that we have incorporated. Hence we have little difficulty in fitting the oscillatory signal from the He ii ionisation zone.

To illustrate the oscillatory signal, Fig. 1 shows a fit to the fourth difference of frequencies for a star with $`M=1M_{}`$ and age $`4.55\times 10^9`$ years. This model is slightly different from the current solar models, since diffusion of helium and heavy elements in the radiative interior is neglected. One can see the two oscillatory components in this figure. The slowly varying component with $`\tau 676`$ s arises from the second helium ionisation zone and the value of $`\tau `$ is approximately equal to the acoustic depth of this ionisation layer. The more rapidly oscillating component ($`\tau 2176`$ s) arises from the base of the convection zone. The amplitudes of all oscillatory components decrease with frequency. At still higher frequencies where the modes are not trapped below the solar surface the effective value of $`\tau `$ for the fits may also change. We do not include such modes as they are not likely to be observed and in any case their interpretation is difficult. Thus for all stellar models we generally restrict ourselves to frequencies less than the acoustic cutoff frequency $`\omega _{ac}\frac{c}{2H_p}`$ where $`c`$ is the sound speed and $`H_p`$ is the pressure scale height near the stellar surface. At low frequencies the fourth difference generally shows steep variations, probably because the low radial order modes do not conform to the asymptotic relations used in obtaining the oscillatory form. Thus we neglect these modes too while studying the oscillatory signal. The frequency range studied for different stars varies with mass and age and in order to compare the respective amplitudes of the oscillatory signal we scale the frequencies by $`\sqrt{M/R^3}`$, where $`M`$ is the stellar mass and $`R`$ the radius in terms of solar values. Thus we compare the amplitudes of oscillatory signal at frequency $`\nu _0\sqrt{\overline{\rho }/\overline{\rho }_{}}`$, where $`\nu _0`$ is a fixed value chosen for the Sun and $`\overline{\rho }`$ is the mean density.

## 3 The stellar models

In order to study the variation in the oscillatory signal for different stars we construct models for stars with different masses at various ages using the stellar evolution code CESAM (Morel 1997). For simplicity we have neglected diffusion of helium and heavy elements in stellar interior. We use the EFF (Eggleton et al. 1973) equation of state with Coulomb correction (Christensen-Dalsgaard & Däppen 1992); the OPAL (Iglesias & Rogers 1996) opacity tables; nuclear reaction rates from the NACRE compilation (Angulo et al. 1999; Morel 1999); mixing length of $`1.8H_p`$, where $`H_p`$ is the pressure scale height.

We have considered stellar models with $`M/M_{}=0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4`$ and $`5`$ at different ages spanning the main sequence phase of evolution. The ZAMS chemical composition of all the models were taken to be solar-like ($`X=0.708,Z=0.019`$). The characteristics of some of the representative models are summarised in Table 1. Here, $`X_c`$ is the hydrogen abundance at the centre, $`\nu _{ac}=c/(4\pi H_p)`$ is the acoustic cutoff frequency at the surface, $`\tau _0`$ is the acoustic radius of the star or the sound travel time from the centre to the surface. The radial distance as well as the acoustic depth of the outer edge of the convective core, base of the convective envelope and the He ii ionisation zone are also given in the table.

For each of these models we calculate the frequencies of p-modes for $`\mathrm{}=0,1,2,3`$ and use them to study the oscillatory signal using the second or fourth difference. The oscillatory signal also depends on the extent of overshoot below the convection zone (Monteiro et al. 1994) as well as the treatment for calculating the diffusion of helium and heavy elements in the radiative regions (Basu & Antia 1994). From the solar studies it has been known that the overshoot below the base of the convection zone is quite small (Basu 1997) and further the amplitude of the observed signal is consistent with that in a solar model without overshoot and with very little composition gradient at the base of the convection zone. The absence of composition gradient is probably due to some mixing process operating in that region (Richard et al. 1996; Brun et al. 1999). Thus in this work we neglect both diffusion and overshoot. Of course, it is possible that in some stars the overshoot may be larger, but that will only increase the amplitude of the oscillatory signal. Similarly, a gradient in helium abundance below the outer convection zone also tends to enhance the oscillatory signal. Thus by neglecting overshoot and diffusion we will get a lower limit on the amplitude of the oscillatory signal. For massive stars where the outer convection zone becomes shallow, it is quite likely that extent of overshoot as well as composition gradient may be larger than that in the Sun, because the efficiency of convection would be lower (Saikia et al. 2000).

Apart from this, for most of these models we also calculate the splitting coefficients assuming a model tachocline profile given by Eq. (3). We use $`\delta \mathrm{\Omega }=20`$ nHz (which is the typical value for the Sun) and $`w=0.001R`$ while $`r_d`$ is taken to be the base of the outer convection zone. There is no established theory to suggest that the tachocline has to coincide with the base of the convection zone, but it is quite likely to be true. In some cases, we study the variation in the oscillatory signal with the depth of the tachocline. The width is chosen to be somewhat small as in massive stars the tachocline may be located in regions where the scale heights are smaller than that for the solar case. The effect of the width of the tachocline on the oscillatory signal has been discussed by Mazumdar & Antia (2001). For stars with $`M1.5M_{}`$ the outer convection zone is too shallow and even a width of $`0.001R`$ may be too large. But such stars are unlikely to have tachoclines and we do not consider this signal in those cases. The variation of amplitude in the splitting coefficients with stellar mass and age is similar to that in the frequencies. We still consider the two separately as the splitting coefficients provide a convenient set to study the variation in oscillatory signal with depth, which in this case can be easily controlled.

## 4 Signal from the base of outer convection zone

In order to study the oscillatory signature in the frequencies we calculate the frequencies for each stellar model and fit the fourth differences of frequencies to an oscillatory signal described in Eq. (4), with two oscillatory components of the same form. The results are shown in Table 2. This table also includes the averaged frequency separations $`\mathrm{\Delta }\nu _0=\nu _{n+1,0}\nu _{n,0}`$, $`D_0=(\nu _{n,0}\nu _{n1,2})/6`$, and $`d_{1/2}=(3/2)(\nu _{n,0}2\nu _{n,1}+\nu _{n+1,0})`$. The fitted $`\tau `$ for the two components are close to the actual acoustic depths of the He ii ionisation zone and the base of the convection zone. Thus it is clear that the depths of these layers can be estimated from the measured frequencies if the accuracy is sufficient to fit the oscillatory signal.

It can be easily shown that taking the second difference of frequencies will enhance the oscillatory signal by a factor of $`4\mathrm{sin}^2(2\pi \tau \mathrm{\Delta }\nu )`$, where $`\mathrm{\Delta }\nu `$ is the frequency separation $`(\nu _{n+1,\mathrm{}}\nu _{n,\mathrm{}})`$. For $`\tau \mathrm{\Delta }\nu <1/12`$ or $`\tau \mathrm{\Delta }\nu >5/12`$ this factor will be smaller than unity and the amplitude of the oscillatory term will reduce when higher order differences are taken. Thus in such cases it may be better to use the second difference or to fit the oscillatory term to the frequencies themselves. In the latter case one will need some model to fit the smooth part also. It can be easily seen that if the smooth part is assumed to be linear in $`n`$, then it will vanish when second difference is taken. Thus we can fit the second difference without worrying about the smooth part. Similarly, if the smooth part is assumed to be a cubic polynomial in $`n`$, it will vanish if the fourth difference is taken. For more general smooth functions, e.g., splines, the smooth part will not necessarily vanish when differences are taken but it will generally become smaller. However, the errors also increase as higher order differences are computed. Thus, fourth differences may be used only if $`\mathrm{sin}^2(2\pi \tau \mathrm{\Delta }\nu )>1/\sqrt{2}`$ (Basu et al. 1994).

If $`\tau _0`$ is the acoustic radius or the sound travel time from the centre to the surface of the star, then $`\mathrm{\Delta }\nu 1/(2\tau _0)`$. Thus the amplification factor for second difference is $`4\mathrm{sin}^2(\pi \tau /\tau _0)`$, which is maximum when $`\tau \tau _0/2`$ or when the layer is halfway between the centre and the surface in terms of sound travel time. On the other hand, if the layer is close to either the surface or the centre, the amplification factor will be small. Thus for massive stars where the outer convection zone is shallow and the convective core is also small in terms of acoustic radius, the fourth difference will reduce the amplitude significantly. In such cases it may be better to use alternate strategy to fit oscillatory signal. In this work since we are using exact frequencies from stellar models the errors are not amplified and hence as long as the amplitude does not decrease significantly when differences are taken we can use the differences to conveniently suppress the smooth part and isolate the oscillatory signal. Even without errors we find that when the outer convection zone becomes too shallow the oscillatory signal reduces significantly when fourth differences are taken. In such cases we have instead fitted the oscillatory signal in second differences. Thus in order to facilitate comparison and remove the variation in amplitude arising from variation of the amplifying factor, we reduce all amplitudes to those in frequencies (or splitting coefficients) by accounting for the amplification factor in taking second or fourth differences. Thus the amplitudes given in Table 2 or those shown in Figs. 2, 3 refer to the frequencies or splitting coefficients, irrespective of whether second or fourth differences were used in the fits. Another advantage of scaling the amplitudes to remove the amplification factor is that variation in the resulting amplitudes is significantly reduced.

Theoretically, there is no reason to believe that taking fourth differences of frequencies will remove the smooth component. But in practice, using a wide range of stellar models, we find that taking the fourth differences essentially removes the smooth component of the frequency in all cases, thus making it much simpler to fit the resulting oscillatory signal. While fitting the second difference or the frequencies directly, it is essential to account for the smooth part, which necessarily involves some assumptions about how the smooth component is represented. The resulting fit will also depend to some extent on these assumptions. In general, the smooth part also depends on $`\mathrm{}`$ which makes it more difficult to remove. For more rapidly varying signal with $`\tau \tau _0/2`$ like that arising from the base of the solar convection zone, it may be relatively easy to characterise the smooth part as has been done by Monteiro et al. (1994). But when the ‘wavelength’ of oscillation is comparable to the frequency range in which observations are available, it is not easy to distinguish between the smooth and oscillatory components. Taking the differences also tends to reduce the amplitude in such cases making it more difficult to fit the oscillatory signal. We believe that this difficulty is more fundamental, as it arises from the fact that there is little difference between the smooth component and an oscillatory component with large wavelength. This problem will manifest itself in some form or other irrespective of whether we fit the differences or the frequencies directly. If taking the fourth differences wipes out the oscillatory signal, then it would be possible to approximate it well by a cubic polynomial, which can define the smooth part.

Fig. 2 shows the variation in amplitude of the fitted signal in the frequencies with stellar age for a few stellar masses. It can be seen that at masses $`<1.5M_{}`$ there is not much variation in amplitude with age and further even the variation with stellar mass is also quite small in the signal from the base of the convective envelope. The maximum difference being about a factor of two between all models included in the figure. The oscillatory signal from the He ii ionisation zone has in general larger amplitude and the amplitude increases with stellar mass. For $`M>M_{}`$ the amplitude of signal from He ii ionisation zone decreases with ages, while at lower mass it is essentially independent of age. For more massive stars ($`M1.5M_{}`$) the acoustic depth of both the He ii ionisation zone and the base of convective envelope are comparable and it becomes difficult to fit the two signals separately. In addition, because the resulting $`\tau `$ is quite small, taking the differences tend to reduce the amplitudes and it becomes more difficult to fit the signal even in the second differences. In some cases, when the frequencies of two components are close the resulting pattern shows beats in the limited frequency range that is likely to be available. The apparent amplitude of the signal depends on the relative phase of the two components and it is difficult to obtain any reliable fit to both the components. The characteristics of the signal arising from the base of the convective envelope may also be different in these stars as the outer convection zone splits into two parts. Consequently, we have not shown the results for these stars.

Apart from the frequencies we also calculate the splitting coefficients for each of the stellar models with a model rotation profile given by Eq. (3). This rotation profile assumes a tachocline at the base of the outer convection zone. The fourth difference of the calculated splitting coefficient $`c_3(n,\mathrm{})`$ for $`\mathrm{}=2,3`$ is again fitted to an oscillatory signature of the form given in Eq. (4) to calculate the characteristics of the tachocline. These results are also included in Table 2 and it is clear that the fitted $`\tau `$ is close to the actual depth of the tachocline in the model profile. The variation of amplitude with stellar mass and age is shown in Fig. 3. In this case, the variation with stellar mass and age is very small.

The typical amplitudes of the oscillatory signal in the frequencies due to the base of outer convection zone is $`0.05\mu `$Hz, while that due to He ii ionisation zone is $`0.5\mu `$Hz. We have not added any errors in the model frequencies and hence there is little difficulty in identifying this signal for low mass ($`M<1.5M_{}`$) stars. Actual observations will naturally have some random error associated with each frequency. From simulations with artificial data sets it turns out that it is possible to detect the oscillatory signal as long as the errors in frequencies are less than or of the order of the amplitude of oscillations. Thus we need an accuracy of about $`0.05\mu `$Hz to detect this signal from the base of convective envelope and measure its characteristics. For $`\tau \tau _0/2`$ the signal to noise ratio will improve by a factor of two when fourth differences are taken and an accuracy of about $`0.1\mu `$Hz may be sufficient to measure the characteristics of the signal. Such an accuracy may be possible with the forthcoming space missions like COROT (Baglin et al. 1998), MOST (Matthews 1998) and MONS (Kjeldsen & Bedding 1998). The oscillatory signal due to the He ii ionisation zone has an order of magnitude larger amplitude for solar mass stars and it would be easier to detect this signal. This signal may be useful in measuring the helium abundance in stellar envelopes. The oscillatory signal in splitting coefficients has much smaller amplitude and unless the stars are rotating much faster and have a much larger $`\delta \mathrm{\Omega }`$ it may not be possible to detect such a signal.

## 5 Signal from the edge of convective core

### 5.1 Signal in frequencies of high mass stars

For stars with $`M1.2M_{}`$, apart from the outer convection zone, a part of the inner core is also convective and we would expect an additional oscillatory component in the frequencies arising from discontinuity in derivatives of sound speed at the outer boundary of the convective core. In this case, $`\tau \tau _0`$ and as mentioned in the previous section, taking fourth difference will reduce the amplitude significantly. In order to understand the oscillatory signal, we consider a stellar model for $`M=2M_{}`$ at an age of $`5\times 10^8`$ years. This model has a convective core with acoustic radius $`0.03\tau _0`$ and taking the fourth difference will reduce the amplitude by a factor of about $`1000`$. The outer convection zone in this model has an acoustic depth of about $`0.08\tau _0`$ and that oscillatory signal will also reduce by a factor of $`20`$. Fig. 4 shows the second and fourth differences of frequencies for this model. It is clear from the figure that the oscillatory signal arising from the outer boundary of the convective core is not seen in either the second or the fourth differences. These differences do show the signal from the base of the outer convective zone or the He ii ionisation region, while the oscillations that are seen in the fourth differences have completely different $`\tau `$. Fitting an oscillatory signal in a restricted frequency range shows that this signal is arising from $`\tau 5300`$ s, which corresponds to $`r/R0.6`$, where we do not expect any boundary of convective region.

In order to understand this signal we show in Fig. 5, the quantity $`W(r)=(1/g)\mathrm{d}c^2/\mathrm{d}r`$ for a few different stellar models. It can be seen that while low mass stars do not have any sharp variation in this scaled derivative of the sound speed in the radiative interior, the high mass stars all show a prominent dip in $`W(r)`$ in the radiative interior. This dip is quite small compared to the sharp rise near the base of the outer convection zone in near-solar mass stars. But for high mass stars where the outer convection zone is very shallow, this dip becomes the dominant feature in the first derivative of the sound speed and hence contributes to the oscillatory signature. Moreover, calculating the fourth difference increases the amplitude of this signal by more than an order of magnitude, while amplitudes of other oscillatory terms coming from edges of convective layers are significantly reduced. Thus in the fourth difference the prominent oscillations are due to this dip in $`W(r)`$, while in the second difference the oscillations arising from the base of the outer convection zone are seen. But the expected signal from the outer edge of the convective core is not seen in either the second or the fourth differences. Since there is no discontinuity in the derivatives of the sound speed in the radiative region the oscillatory signal does not have the same form as that expected from a discontinuity, but for low frequency modes, where the radial wavelength would be larger than the size of the dip the signal is seen more clearly. At high frequencies the radial wavelength may be comparable to or smaller than the size of the dip in $`W(r)`$ and the oscillatory signal is not clear. Further, just like the He ii ionisation zone, the amplitude of this oscillatory signal should also be modulated by an oscillatory factor, which is not included in our fits.

### 5.2 Signal in splitting coefficients

In order to understand our failure to see any oscillatory signal from the convective core in the frequencies, we attempt to study the oscillatory signal due to the tachocline in the splitting coefficients. We calculate the splitting coefficients assuming the tachocline to have different depths to study the variation in amplitude with depth and the results are shown in Table 3. We have considered a stellar model with mass $`2M_{}`$ and age $`0.5`$ Gyr for the purpose of illustration. The amplitude appears to be maximum when the tachocline is located near the acoustic midpoint of the star (in this particular case, near $`r/R=0.8`$) and falls off as it shifts to either end. It can be easily seen that this variation is more or less fully explained by the amplification factor while taking the differences. Thus it appears that the amplitude of the signal in the splitting coefficients is more or less independent of the depth and the apparent variations are coming from calculating the differences. It is clear that in order to detect the oscillatory signal from the outer edge of the convective core, we should look at the splitting coefficients directly, rather than their differences.

Fig. 6 shows the second difference of splitting coefficients for $`\mathrm{}=2`$ modes arising from the tachocline at different depths in the stellar model mentioned above. It can be seen that the wavelength of oscillations reduces starting from $`r_d=0.95R`$ as we go deeper as expected. However, after $`r_d=0.8R`$ the wavelength appears to increase again contrary to what we expect. This is similar to the phenomenon of aliasing in Fourier transform. For simplicity if we assume that the frequencies of $`\mathrm{}=2`$ modes are uniformly spaced then the oscillatory term can be written as

$`c_3(n,\mathrm{})`$ $`=`$ $`A\mathrm{sin}[4\pi \tau (\nu _0+n\mathrm{\Delta }\nu _{\mathrm{}})+\varphi ]`$ (5)
$`=`$ $`A\mathrm{sin}[4\pi \tau \nu _0+\varphi +2\pi n\tau /\tau _0],`$
where $`\mathrm{\Delta }\nu _{\mathrm{}}`$ is the spacing between successive modes for given $`\mathrm{}`$. Here we have used the fact that for low $`\mathrm{}`$, $`\mathrm{\Delta }\nu _{\mathrm{}}1/(2\tau _0)`$. Now if $`\tau =\tau _0`$ it is clear that for all $`n`$ the oscillatory term will be the same and $`c_3(n,\mathrm{})`$ will not show any oscillations. If we take $`\tau ^{}=\tau _0\tau `$, then it is clear that the oscillatory term will be the same apart from a change in the phase for all modes. Thus from the oscillatory signal for an individual $`\mathrm{}`$ it is not possible to distinguish between these two values of $`\tau `$. Fig. 6 shows the actual fit obtained using two different $`\tau `$ for some depths. This is exactly similar to aliasing in discrete Fourier transform. Alternatively, it means that we can measure the acoustic distance from either boundary. If we consider sampling at a uniform spacing of $`\mathrm{\Delta }\nu 1/(2\tau _0)`$, then the corresponding Nyquist frequency will be $`1/(2\mathrm{\Delta }\nu )\tau _0`$. Thus ‘frequencies’ higher than $`\tau _0`$ will be aliased to lower frequencies. Because of the factor of $`2`$ inside the sine function in Eq. (4) this actually corresponds to an acoustic depth of $`\tau _0/2`$. Thus this is identical to aliasing in discrete Fourier transform.

It may be noted that the lowest panel in Fig. 6 for $`r_d=0.1R`$ is shown separately, as in this case instead of the second differences (which have very small amplitude), we show the splitting coefficients themselves. Although, the amplitude of variation in this is comparable to amplitudes at other depths, it is impossible to distinguish this signal from a smooth trend. At other depths we have shown the second difference as the presence of smooth trend in splitting coefficients does not show the oscillatory signal clearly. It is clear from this panel that both $`\mathrm{}=2`$ and $`3`$ modes separately fall on a smooth curve and only if they are combined together the oscillatory signal from the core can be identified. It may be noted that the form of rotation rate used to calculate these splitting coefficients does not have any smooth component and hence the smooth component in resulting splitting coefficients is quite small, thus enabling us to see the signal in some form. In real stars there will also be some smooth component which will be difficult to isolate from the oscillatory part when the apparent $`\tau `$ is quite small.

In principle, the ambiguity in $`\tau `$ can be resolved by including multiple $`\mathrm{}`$ values. The two panels at the lower end of Fig. 6 shows the splitting coefficients for $`r_d=0.3R`$ and $`0.1R`$ for $`\mathrm{}=2,3`$. It is clear that although modes corresponding to each $`\mathrm{}`$ show a variation with large wavelength, the two do not agree with each other in phase and if we attempt a fit using Eq. (4), the fit will not be good when $`\tau `$ corresponds to the smaller value implied by individual $`\mathrm{}`$. Only an oscillatory form with larger value of $`\tau `$ can fit all the modes. In terms of discrete Fourier transform this can be understood as follows: because of additional $`\mathrm{}`$ the sampling interval is reduced thus increasing the Nyquist frequency, enabling us to determine higher frequencies. However, the problem with fitting any form is that in all cases there will be some smooth component in the splitting coefficient arising from smooth variations in the rotation rate and it will be difficult to separate out this component which has to be done for each $`\mathrm{}`$ separately, as the smooth component will depend on $`\mathrm{}`$ and frequency. The problem arises because in this case the oscillatory component also appears to be smooth due to aliasing and it will be difficult to isolate it from the smooth component as can be seen by the panel for $`r_d=0.1R`$ in Fig. 6. This separation will be more difficult when errors are added to the splitting coefficients, as will be the case in real observed values. Thus it is difficult to detect any possible oscillatory signature arising from convective cores in high mass stars.

### 5.3 Signal from the convective core

Detecting the oscillatory signal in the frequency from the convective core is even more difficult as the frequencies have dominant oscillatory contributions arising from other variations in sound speed derivatives. Thus the signal due to the outer convection zone will have similar $`\tau `$ to that from convective core (after aliasing) and it will be difficult to isolate the core contribution from the surface contribution. Further, there is a smooth variation with $`n`$ for the frequencies, which is several orders of magnitude larger than the oscillatory part and it will be more difficult to remove this smooth variation. Of all the models tried we could detect this signal in only $`1.75M_{}`$ star with an age around $`0.9`$ Gyr (Fig. 7). Interestingly even with exact frequency data this signal is not clearly seen at all ages. We believe the main reason for this is that at this particular age the amplitude of signal arising from the outer convection zone and He ii ionisation zone combine to yield a small amplitude, possibly because of beating between the two frequencies which are reasonably close. Because of the small amplitude, even smaller oscillatory signal due to the convective core shows up with amplitude of order of $`0.01\mu `$Hz, which would otherwise have been missed. This signal will translate to an amplitude of $`0.3\mu `$Hz in frequencies, which can be compared with frequency separation of $`50\mu `$Hz. Although, the amplitude of this oscillatory signal is fairly large in the frequencies themselves it gets reduced when differences are taken. If we do not take differences we need to remove the smooth part independently, which is also difficult as explained earlier. We also tried to use the first difference of frequency instead of second difference, so that the amplitude of oscillatory signal from the convective core would be larger. But even in that case the smooth part is substantial and it is not easy to remove it in most cases, thus making it difficult to identify the oscillatory signal.

In an attempt to detect the oscillatory signature we consider even higher mass stars where the convective core is comparatively larger. Thus we consider a stellar model with $`M=5M_{}`$ and age $`5\times 10^7`$ years. In this case the convective core has a radius of $`0.04\tau _0`$ and Fig. 7 shows the second differences of the frequencies. It appears that at a somewhat high frequency range some oscillatory signal which is similar to what one would expect from the convective core is present. However, the acoustic cutoff frequency for this model is only $`0.72`$ mHz and hence the oscillatory pattern is seen only in high frequency modes which are not trapped in the interior. These modes are unlikely to be seen in distant stars. Even if they are seen, it will be difficult to determine their frequencies accurately, because of large widths in the power spectrum that may be expected. Similar signal at frequencies above the acoustic cutoff has been seen in some other stellar models also. The amplitude of dominant oscillations due to the He ii ionisation zone falls off rapidly with frequency and that probably allows the smaller signal arising from the convective core to be seen at high frequencies.

Since we have failed to detect the oscillatory signal from convective core, in the next section we propose an alternative technique which looks at the variations in the smooth part to detect convective cores in stars and possibly measure its size.

## 6 Signature of convective core in frequency separation

From the asymptotic theory of stellar oscillations it is well known that the frequency separations

$`\mathrm{\Delta }\nu _{\mathrm{}}`$ $`=`$ $`\nu _{n+1,\mathrm{}}\nu _{n,\mathrm{}},`$
$`D_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{\nu _{n,\mathrm{}}\nu _{n1,\mathrm{}+2}}{4\mathrm{}+6}},`$ (6)
$`d_{1/2}`$ $`=`$ $`{\displaystyle \frac{3}{2}}(\nu _{n,0}2\nu _{n,1}+\nu _{n+1,0}),`$
are related to variation in the sound speed across the stellar radius including the convective core. Some of these frequency separations are quite sensitive to conditions in stellar core. These frequency separations can be easily measured for stars using only the low degree modes. Here we consider only the frequency separations for $`\mathrm{}=0`$ averaged over a frequency range scaled appropriately for each model to correspond to a range of $`3`$$`4`$ mHz for the Sun.

Christensen-Dalsgaard (1993) has suggested that these frequency separations can be related to the stellar age and mass from the evolutionary tracks in the ($`\mathrm{\Delta }\nu _0`$,$`D_0`$) plane, which is also known as the C-D diagram. Since the frequencies of stellar oscillations are expected to scale as $`\sqrt{M/R^3}`$ it will be more instructive to scale these frequency differences by this factor to isolate other effects in the frequency variation. Fig. 9 shows the scaled frequency separations plotted against each other for various stellar models. This is similar to the figure shown by Christensen-Dalsgaard (1998). It can be immediately seen that after scaling the variation in the frequency separation is significantly reduced and moreover, all models for $`M1.25M_{}`$ appear to have the scaled large frequency separation $`\mathrm{\Delta }\nu _0/\sqrt{\overline{\rho }/\overline{\rho }_{}}137\mu `$Hz, close to the solar value. These are the models that do not have any significant convective core. All models with convective core appear to have significantly smaller scaled frequency separations. Thus it appears that if the stellar mass and radius are known then by looking at the scaled frequency separation we can infer the presence of a convective core (cf., Thompson 2000). If the scaled large frequency separation is less than about $`134\mu `$Hz, then the star is likely to have a convective core.

Apart from detecting the presence of convective core we would also like to measure its size. Indeed, the asymptotic theory of stellar oscillations (e.g., Roxburgh & Vorontsov 2001) may be able to provide a relation connecting the frequency separation to the size of the convective core. In this work, using the frequencies computed from different stellar models we attempt to check if the deviation from the solar value of the large frequency separation can be correlated to the size of the convective core. It might also be possible to derive some combination of the large and small frequency separations which would vary in tandem with the size of the convective core. In Fig. 10 we have plotted the difference of the scaled large frequency separation, from the solar value of $`137\mu `$Hz, against the radius of the convective core for stellar models of different mass and age. Both the axes have been scaled to bring the curves corresponding to different masses close to each other. We find that there is an approximately linear relation between these two quantities for a wide range of masses. Although, all the curves do not collapse into one single relationship, this diagram might be used to estimate the size of the convective core, provided the mass and radius of the star are independently known. Only the curves for $`1.25M_{}`$ and $`1.5M_{}`$ stars appear to fall outside the band.

The small frequency separation, $`D_0`$, is found to bear a correlation with the radial size of the convective core. This relation becomes clear with appropriate scaling of $`D_0`$ as well as $`r_{\mathrm{core}}`$, as shown in Fig. 11(a). One problem with using $`D_0`$ to measure the convective core size is that the range of values covered by a solar mass star during its evolution has a large overlap with the value for stars with convective core. Thus by looking at the value of $`D_0`$ alone it will be difficult to determine the convective core size.

Basu et al. (2001) have pointed out that the behaviour of the small separation, $`d_{1/2}`$, is somewhat different from the other frequency separations. We find that this quantity does not have a clear relationship with the size of the core. Nevertheless, as shown in Fig. 11(b), its value derived for models with a convective core (after some scaling) is conspicuously different from that of models without convective core. However, for both the small separations, we are unable to collapse all the curves together as would have been desirable for them to be used as definitive diagnostics for the size of a convective core. However, they can be used as complementary diagnostics for properties of convective core.

An advantage of using the frequency separation as opposed to the oscillatory signal discussed in the previous section is that comparatively lower accuracy is required in this case. Thus an accuracy of order of $`1\mu `$Hz may be sufficient to get an estimate of convective core size using the frequency separations. On the other hand, extracting the oscillatory signal will require at least an order of magnitude higher accuracy in measured frequencies. Even with very high accuracy it will be difficult to isolate the oscillatory signal from convective core from other signals. However, the use of frequency separation requires the knowledge of stellar mass and radius for scaling the frequencies. This may not be a problem for nearby stars which are likely to be the prime targets of asteroseismic missions.

## 7 Conclusions

In this work we have tried to study the oscillatory signal in the frequency arising from discontinuities in derivatives of sound speed inside a star. This oscillatory signal can be conveniently used to study the location of the base of the outer convection zone or the He ii ionisation zone as has been suggested by Monteiro et al. (2000). Our results also support their conclusions. We have studied this signal in a wide range of stellar models covering a range of stellar masses and ages. One of the reasons for this study was to possibly identify promising stars where such signal may be large enough to be detected. However, from our study it appears that there is only a small variation in the amplitude of the oscillatory signal arising from the base of the convection zone between stars of different masses and ages on the main sequence. The amplitude increases with decreasing mass or age. The amplitude of the oscillatory signal arising from the He ii ionisation zone appears to increase with stellar mass and can possibly be used to measure helium abundance in stellar envelopes. Typical amplitude of oscillations from the base of the outer convection zone is around $`0.05\mu `$Hz and from the simulations that we have carried out it appears that comparable accuracy will be required to extract this signal and measure the characteristics of the convection zone boundary. Although, in this study, we have not included any overshoot below the convection zone, it is well known that the amplitude of this signal increases with overshoot and it will be possible to use this signal to measure the extent of overshoot below the base of the convection zone (Monteiro et al. 2000). For the case of the Sun, this has already been done, though with the help of higher degree modes whose frequencies are known more accurately. The oscillatory signal arising from the He ii ionisation zone is an order of magnitude larger and it should be possible to detect this signal in observed frequencies.

In higher mass stars we find an oscillatory signal arising from the dip in $`W(r)`$ in the radiative interior, which may also provide some diagnostic for stellar models. From Fig. 5 it can be seen that the location of the dip depends on stellar mass. Since this location can be determined through $`\tau `$ in oscillatory signal it can provide some constraints on stellar models. It can be shown that in radiative regions $`W(r)`$ is also affected by opacity, which in turn will depend on heavy element abundance, $`Z`$. Thus the depth of the dip also depends on $`Z`$. However, the amplitude of oscillations arising from this dip is about $`0.05\mu `$Hz in the fourth differences and hence we will need accuracy of order of $`0.01\mu `$Hz in measured frequencies to detect this signal. This level of accuracy may not be achieved in immediate future, but with further improvements and longer observations it may be possible to detect this signal.

Most of the earlier studies have concentrated on the oscillatory signal from the outer convection zone. It is clear that there would be a discontinuity in the derivative of sound speed at the outer edge of the convective core in a massive star. This discontinuity should also give rise to a similar oscillatory signal. Although, Monteiro et al. (1998) have suggested that this signal can be used to study characteristics of convective core, it is not clear if this signal has been seen in frequencies of a stellar model. Thus in this work we attempt to study the oscillatory signature arising from a convective core. Even using exact frequencies from stellar models we find it difficult to detect this signal. We have tried to identify the reasons for this failure by studying a similar signal in the splitting coefficients arising from a possible tachocline at different depths. We find that the oscillations display a phenomenon similar to aliasing in discrete Fourier transform, because of which using modes from a single value of $`\mathrm{}`$, it is not possible to distinguish between oscillations arising from discontinuities at acoustic depths of $`\tau `$ and $`\tau _0\tau `$, where $`\tau _0`$ is the total acoustic radius of the star. Due to this aliasing the signal from the convective core would appear as a smooth variation in frequencies from any individual $`\mathrm{}`$. Though, in principle, by combining data for $`\mathrm{}=0,1,2,3`$ we should be able to characterise this signal, but in practice we have to remove the smooth part for each $`\mathrm{}`$ before the oscillatory part can be studied. It is not easy to separate this signal from the convective core, thus making it difficult to study its characteristics.

From our attempts to find the oscillatory signal from the convective core, we could find this signature in a few stellar models at frequencies above the acoustic cutoff frequencies. Such modes are not likely to be observed in stars, but from these fits it appears that the amplitude of the oscillatory signal in frequencies is actually quite large, around $`0.3\mu `$Hz, but because of our inability to separate it from the smooth component it is difficult to detect this signal. Of all the models that we tried we could fit this signal only in a $`1.75M_{}`$ star around the age of $`0.9`$ Gyr at frequencies less than the acoustic cutoff frequency. This is most probably because the dominant oscillatory signal due to outer convection zone and He ii ionisation zone combine to yield a small amplitude in this stellar model, because of beating. Thus we conclude that an oscillatory signal of the form studied in this work does not provide a viable means of studying convective cores in massive stars.

Having failed to detect the oscillatory signature from the convective core in frequencies of stellar oscillations, we attempt to check if the smooth part of the frequency can be used to study the properties of the core. The frequency separations $`\mathrm{\Delta }\nu _{\mathrm{}}`$, $`D_{\mathrm{}}`$ and $`d_{1/2}`$ defined by Eq. (6) are known to be sensitive to stellar structure including the core (e.g., Christensen-Dalsgaard 1993; Basu et al. 2001). Since the frequencies of stellar oscillations generally scale with $`\sqrt{M/R^3}`$ it is instructive to consider the scaled frequency separations. It turns out that the scaled large frequency separation is generally around $`137\mu `$Hz for stars without convective core, while for massive stars with convective cores it tends to be significantly smaller. Thus this scaled frequency separation can be used to detect the presence of a convective core. After suitable scaling the departure of this frequency separation from the mean value for stars without a convective core can be roughly correlated to the size of the convective core as can be seen from Fig. 10. It can be seen that apart from intermediate masses of $`1.25M_{}`$ and $`1.5M_{}`$, for higher masses all the evolutionary tracks lie in a band and from the measure of the frequency separation together with the scaling factor it should be possible to estimate an approximate size of the convective core. The small frequency separations $`D_0`$ and $`d_{1/2}`$ are also sensitive to properties of the core, though we could not find a similar tight relation with the core size. We believe that by combining all the information in the smooth component of the frequency for only low degree modes it should be possible to study some properties of the convective core, including its size. An advantage of using these frequency separations is that comparatively lower accuracy of $`1\mu `$Hz in frequency would be sufficient, while we would need much higher accuracy of about $`0.05\mu `$Hz to study the oscillatory signal. Even the oscillatory signal from the base of the solar convection zone has been mostly studied using intermediate degree modes which have higher accuracy. In this work we have considered only the mean frequency separations which naturally discards some information, which is contained in the individual frequency separations for each value of $`n`$. In particular, the variation of frequency separation with $`n`$ will also give some useful information about the convective core. Instead of frequency separation inversion techniques using only low degree modes have also been tried to study stellar cores (Basu et al. 2001).

In the foregoing discussion it is assumed that the scaling factor $`\sqrt{M/R^3}`$ is known for stars. This would require knowledge of stellar mass and radius. In principle, the stellar mass and age can be determined from the C-D diagram (Christensen-Dalsgaard 1993), while if the distance to the star is known, its luminosity and hence the radius can be determined. From Fig. 9 it can be seen that the scaled large frequency separation is not very sensitive to stellar mass or age. The variation over the entire range of stellar models considered by us is less than 10%. Hence in absence of any other information the measured large frequency separation can also be used to measure the scale factor. But in that case the resulting separation cannot be used to measure the size of convective core.

## Acknowledgements

We thank P. Morel for providing the CESAM code for stellar evolution and for the help in its usage, and S. M. Chitre for useful discussions. We also thank the Referee, M. J. P. F. G. Monteiro for useful suggestions.


format: text

text:
# Quantum Clustering

## 1 Introduction

Clustering of data is a well-known problem of pattern recognition, covered in textbooks such as . The problem we are looking at is defining clusters of data solely by the proximity of data points to one another. This problem is one of unsupervised learning, and is in general ill defined. Solutions to such problems can be based on intuition derived from Physics. A good example of the latter is the algorithm by  that is based on associating points with Potts-spins and formulating an appropriate model of statistical mechanics. We propose an alternative that is also based on physical intuition, this one being derived from quantum mechanics.

As an introduction to our approach we start with the scale-space algorithm by  who uses a Parzen-window estimator of the probability distribution leading to the data at hand. The estimator is constructed by associating a Gaussian with each of the $`N`$ data points in a Euclidean space of dimension $`d`$ and summing over all of them. This can be represented, up to an overall normalization by

$$\psi (𝐱)=\underset{i}{}e^{\frac{(𝐱𝐱_i)^2}{2\sigma ^2}}$$
(1)
where $`𝐱_i`$ are the data points. Roberts  views the maxima of this function as determining the locations of cluster centers.

An alternative, and somewhat related method, is Support Vector Clustering (SVC)  that is based on a Hilbert-space analysis. In SVC one defines a transformation from data space to vectors in an abstract Hilbert space. SVC proceeds to search for the minimal sphere surrounding these states in Hilbert space.We will also associate data-points with states in Hilbert space. Such states may be represented by Gaussian wave-functions, in which case $`\psi (𝐱)`$ is the sum of all these states. This is the starting point of our method. We will search for the Schrödinger potential for which $`\psi (𝐱)`$ is a ground-state. The minima of the potential define our cluster centers.

## 2  The Schrödinger Potential

We wish to view $`\psi `$ as an eigenstate of the Schrödinger equation

$$H\psi (\frac{\sigma ^2}{2}^2+V(𝐱))\psi =E\psi .$$
(2)
Here we rescaled $`H`$ and $`V`$ of the conventional quantum mechanical equation to leave only one free parameter, $`\sigma `$. For comparison, the case of a single point at $`𝐱_1`$ corresponds to Eq. 2 with $`V=\frac{1}{2\sigma ^2}(𝐱𝐱_1)^2`$ and $`E=d/2`$, thus coinciding with the ground state of the harmonic oscillator in quantum mechanics.

Given $`\psi `$ for any set of data points we can solve Eq. 2 for $`V`$:

$$V(𝐱)=E+\frac{\frac{\sigma ^2}{2}^2\psi }{\psi }=E\frac{d}{2}+\frac{1}{2\sigma ^2\psi }\underset{i}{}(𝐱𝐱_i)^2e^{\frac{(𝐱𝐱_i)^2}{2\sigma ^2}}$$
(3)
Let us furthermore require that min$`V`$=0. This sets the value of

$$E=\mathrm{min}\frac{\frac{\sigma ^2}{2}^2\psi }{\psi }$$
(4)
and determines $`V(𝐱)`$ uniquely. $`E`$ has to be positive since $`V`$ is a non-negative function. Moreover, since the last term in Eq. 3 is positive definite, it follows that

$$0<E\frac{d}{2}.$$
(5)
We note that $`\psi `$ is positive-definite. Hence, being an eigenfunction of the operator $`H`$ in Eq. (2), its eigenvalue $`E`$ is the lowest eigenvalue of $`H`$, i.e. it describes the ground-state. All higher eigenfunctions have nodes whose number increases as their energy eigenvalues increase.111In quantum mechanics, where one interprets $`|\psi |^2`$ as the probability distribution, all eigenfunctions of $`H`$ have physical meaning. Although this approach could be adopted, we have chosen $`\psi `$ as the probability distribution because of simplicity of algebraic manipulations.

Given a set of points defined within some region of space, we expect $`V(𝐱)`$ to grow quadratically outside this region, and to exhibit one or several local minima within the region. Just as in the harmonic potential of the single point problem, we expect the ground-state wave-function to concentrate around the minima of the potential $`V(𝐱)`$. Therefore we will identify these minima with cluster centers.

As an example we display results for the crab data set taken from Ripley’s book . These data, given in a five-dimensional parameter space, show nice separation of the four classes contained in them when displayed in two dimensions spanned by the 2nd and 3rd principal components (eigenvectors) of the correlation matrix of the data. The information supplied to the clustering algorithm contains only the coordinates of the data-points. We display the correct classification to allow for visual comparison of the clustering method with the data. Starting with $`\sigma =1/\sqrt{2}`$ we see in Fig. 1 that the Parzen probability distribution, or the wave-function $`\psi `$, has only a single maximum. Nonetheless, the potential, displayed in Fig. 2, shows already four minima at the relevant locations. The overlap of the topographic map of the potential with the true classification is quite amazing. The minima are the centers of attraction of the potential, and they are clearly evident although the wave-function does not display local maxima at these points. The fact that $`V(𝐱)=E`$ lies above the range where all valleys merge explains why $`\psi (𝐱)`$ is smoothly distributed over the whole domain.

We ran our method also on the iris data set , which is a standard benchmark obtainable from the UCI repository . The data set contains 150 instances each composed of four measurements of an iris flower. There are three types of flowers, represented by 50 instances each. Clustering of these data in the space of the first two principal components is depicted in Figure 3. The results shown here are for $`\sigma =1/4`$ which allowed for clustering that amounts to misclassification of three instances only, quite a remarkable achievement (see, e.g. ).

## 3  Application of Quantum Clustering

The examples displayed in the previous section show that, if the spatial representation of the data allows for meaningful clustering using geometric information, quantum clustering (QC) will do the job. There remain, however, several technical questions to be answered: What is the preferred choice of $`\sigma `$? How can QC be applied in high dimensions? How does one choose the appropriate space, or metric, in which to perform the analysis? We will confront these issues in the following.

In the crabs-data of we find that as $`\sigma `$ is decreased to $`\frac{1}{2}`$, the previous minima of $`V(𝐱)`$ get deeper and two new minima are formed. However the latter are insignificant, in the sense that they lie at high values (of order $`E`$), as shown in Fig. 4. Thus, if we classify data-points to clusters according to their topographic location on the surface of $`V(𝐱)`$, roughly the same clustering assignment is expected for $`\sigma =\frac{1}{2}`$ as for $`\frac{1}{\sqrt{2}}`$. One important advantage of quantum clustering is that $`E`$ sets the scale on which minima are observed. Thus, we learn from Fig. 2 that the cores of all four clusters can be found at $`V`$ values below $`.4E`$. The same holds for the significant minima of Fig. 4.

By the way, the wave function acquires only one additional maximum at $`\sigma =\frac{1}{2}`$. As $`\sigma `$ is being further decreased, more and more maxima are expected in $`\psi `$ and an ever increasing number of minima (limited by $`N`$) in $`V`$.

The one parameter of our problem, $`\sigma `$, signifies the distance that we probe. Accordingly we expect to find clusters relevant to proximity information of the same order of magnitude. One may therefore vary $`\sigma `$ continuously and look for stability of cluster solutions, or limit oneself to relatively low values of $`\sigma `$ and decide to stop the search once a few clusters are being uncovered.

## 4  Principal Component Metrics

In all our examples data were given in some high-dimensional space and we have analyzed them after defining a projection and a metric, using the PCA approach. The latter defines a metric that is intrinsic to the data, determined by their second order statistics. But even then, several possibilities exist, leading to non-equivalent results.

Principal component decomposition can be applied both to the correlation matrix $`C_{\alpha \beta }=<x_\alpha x_\beta >`$ and to the covariance matrix

$$𝒞_{\alpha \beta }=<(x_\alpha <x>_\alpha )(x_\beta <x>_\beta )>=C_{\alpha \beta }<x>_\alpha <x>_\beta .$$
(6)
In both cases averaging is performed over all data points, and the indices indicate spatial coordinates from 1 to $`d`$. The principal components are the eigenvectors of these matrices. Thus we have two natural bases in which to represent the data. Moreover one often renormalizes the eigenvector projections, dividing them by the square-roots of their eigenvalues. This procedure is known as “whitening”, leading to a renormalized correlation or covariance matrix of unity. This is a scale-free representation, that would naturally lead one to start with $`\sigma =1`$ in the search for (higher-order) structure of the data.

The PCA approach that we have used on our examples was based on whitened correlation matrix projections. This turns out to produce good separation of crab-data in PC2-PC3 and of iris-data in PC1-PC2. Had we used the covariance matrix $`𝒞`$ instead, we would get similar, but slightly worse, separation of crab-data, and a much worse representation of the iris-data. Our examples are meant to convince the reader that once a good metric is found, QC conveys the correct information. Hence we allowed ourselves to search first for the best geometric representation, and then apply QC.

## 5  QC in Higher Dimensions

In the iris problem we obtained excellent clustering results using the first two principal components, whereas in the crabs problem, clustering that depicts correctly the classification necessitates components 2 and 3. However, once this is realized, it does not harm to add the 1st component. This requires working in a 3-dimensional space, spanned by the three leading PCs. Increasing dimensionality means higher computational complexity, often limiting the applicability of a numerical method. Nonetheless, here we can overcome this “curse of dimensionality” by limiting ourselves to evaluating $`V`$ at locations of data-points only. Since we are interested in where the minima lie, and since invariably they lie near data points, no harm is done by this limitation. The results are depicted in Fig. 5. Shown here are $`V/E`$ values as function of the serial number of the data, using the same symbols as in Fig. 2 to allow for visual comparison. Using all data of $`V<0.3E`$ one obtains cluster cores that are well separated in space, corresponding to the four classes that exist in the data. Only 9 of the 129 points that obey $`V<0.3E`$ are misclassified by this procedure. Adding higher PCs, first component 4 and then component 5, leads to deterioration in clustering quality. In particular, lower cutoffs in $`V/E`$, including lower fractions of data, are required to define cluster cores that are well separated in their relevant spaces.

One may locate the cluster centers, and deduce the clustering allocation of the data, by following dynamics of gradient descent into the potential minima. Defining $`𝐲_i(0)=𝐱_i`$ one follows steps of $`𝐲_i(t+\mathrm{\Delta }t)=𝐲_i(t)\eta (t)V(𝐲_i(t))`$, letting the points $`𝐲_i`$ reach an asymptotic fixed value coinciding with a cluster center. More sophisticated minimum search algorithms (see, e.g. chapter 10 in ) can be applied to reach the fixed points faster. The results of a gradient-descent procedure, applied to the 3D analysis of the crabs data shown in Fig. 5, is displayed in Fig. 6. One clearly observes the four clusters, and can compare clustering allocation with the original data classes.

## 6  Distance-based QC Formulation

Gradient descent calls for the calculation of $`V`$ both on the original data-points, as well as on the trajectories they follow. An alternative approach can be to restrict oneself to the original values of $`V`$, as in the example displayed in Figure 5, and follow a hybrid algorithm to be described below. Before turning to such an algorithm let us note that, in this case, we evaluate $`V`$ on a discrete set of points $`V(𝐱_i)=V_i`$. We can then express $`V`$ in terms of the distance matrix $`D_{ij}=|𝐱_i𝐱_j|`$ as

$$V_i=E\frac{d}{2}+\frac{1}{2\sigma ^2}\frac{\underset{j}{}D_{ij}^2e^{\frac{D_{ij}^2}{2\sigma ^2}}}{_je^{\frac{D_{ij}^2}{2\sigma ^2}}}$$
(7)
with $`E`$ chosen appropriately so that min$`V_i`$=0. This kind of formulation is of particular importance if the original information is given in terms of distances between data points rather than their locations in space. In this case we have to proceed with distance information only.

Applying QC we can reach results of the type of Fig. 5 without invoking any explicit spatial distribution of the points in question. One may then analyze the results by choosing a cutoff, e.g. $`V<0.2E`$, such that a fraction (e.g. $`1/3`$) of the data will be included. On this subset select groups of points whose distances from one another are smaller than, e.g., 2$`\sigma `$, thus defining cores of clusters. Then continue with higher values of $`V`$, e.g. $`0.2E<V<0.4E`$, allocating points to previous clusters or forming new cores. Since the choice of distance cutoff in cluster allocation is quite arbitrary, this method cannot be guaranteed to work as well as the gradient-descent approach.

## 7  Generalization

Our method can be easily generalized to allow for different weighting of different points, as in

$$\psi (𝐱)=\underset{i}{}c_ie^{\frac{(𝐱𝐱_i)^2}{2\sigma ^2}}$$
(8)
with $`c_i0`$. This is important if we have some prior information or some other means for emphasizing or deemphasizing the influence of data points. An example of the latter is using QC in conjunction with SVC . SVC has the possibility of labelling points as outliers. This is done by applying quadratic maximization to the Lagrangian

$$W=1\underset{i,j}{}\beta _i\beta _je^{\frac{(𝐱_i𝐱_j)^2}{2\sigma ^2}}$$
(9)
over the space of all $`0\beta _i\frac{1}{pN}`$ subject to the constraint $`_i\beta _i=1`$. The points for which the upper bound of $`\beta _i`$ is reached are labelled as outliers. Their number is regulated by $`p`$, being limited by $`pN`$. Using for the QC analysis a choice of $`c_i=\frac{1}{pN}\beta _i`$ will eliminate the outliers of SVC and emphasize the role of the points expected to lie within the clusters.

## 8  Discussion

QC constructs a potential function $`V(𝐱)`$ on the basis of data-points, using one parameter, $`\sigma `$, that controls the width of the structures that we search for. The advantage of the potential $`V`$ over the scale-space probability distribution is that the minima of the former are better defined (deep and robust) than the maxima of the latter. However, both of these methods put the emphasis on cluster centers, rather than, e.g., cluster boundaries. Since the equipotentials of $`V`$ may take arbitrary shapes, the clusters are not spherical, as in the k-means approach. Nonethelss, spherical clusters appear more naturally than, e.g., ring-shaped or toroidal clusters, even if the data would accomodate them. If some global symmetry is to be expected, e.g. global spherical symmetry, it can be incorporated in the original Schrödinger equation defining the potential function.

QC can be applied in high dimensions by limiting the evaluation of the potential, given as an explicit analytic expression of Gaussian terms, to locations of data points only. Thus the complexity of evaluating $`V_i`$ is of order $`N^2`$ independently of dimensionality. Our algorithm has one free parameter, the scale $`\sigma `$. In all examples we confined ourselves to scales that are of order 1, because we have worked within whitened PCA spaces. If our method is applied to a different data-space, the range of scales to be searched for could be determined by some other prior information.

Since the strength of our algorithm lies in the easy selection of cluster cores, it can be used as a first stage of a hybrid approach employing other techniques after the identification of cluster centers. The fact that we do not have to take care of feeble minima, but consider only robust deep minima, turns the identification of a core into an easy problem. Thus, an approach that drives its rationale from physical intuition in quantum mechanics, can lead to interesting results in the field of pattern classification.

We thank B. Reznik for a helpful discussion.


text:
# Colloidal dipolar interactions in 2D smectic C films

## I Introduction

Colloidal systems consist of a dispersive medium or solvent and a dispersed phase. The size of the dispersed domains or colloids is large compared to atoms but small compared to macroscopic lengths. Thus these systems exhibit new phenomena that are challenging theoretically while they are of practical importance in a variety of applications, from paints and coatings to food and drugs .

The colloidal state is metastable and its life-time is determined (largely) by the effective interactions between colloids. In many systems these interactions are tunable through the control of one or more external parameters. To prevent phase separation due to van der Waals or other attractive interactions, the colloids are coated or charged to induce steric or electrostatic repulsions . Recently, a novel mechanism of colloidal interactions was reported: the elastic distortion of a liquid crystal host in a reverse nematic emulsion , i.e., an emulsion of water droplets in a nematic liquid. This interaction arises from a competition between the aligning properties of the spherical droplet surfaces that favour e.g., a radial (homeotropic) nematic orientation and the bulk elasticity and boundary conditions that favour a uniform nematic. This effect is generic and it arises in a variety of anisotropic liquids for different isotropic inclusions. In addition, elastic mediated interactions were shown to be effective above the bulk nematic-isotropic temperature, for colloids wet by the nematic phase . Finally, under appropriate conditions, inverted nematic emulsions exhibit colloidal self-organization, i.e., new phases where the colloids are ordered over macroscopic distances .

Understanding the phase behaviour of these complex systems begins by identifying the mechanisms and calculating the effective interactions between colloids. This is far from a simple task owing to the non-linear nature of the nematic elasticity and the existence of widely different length scales: the bulk correlation length that sets the scale for variations of the microscopic degrees of freedom, such as nematic order, biaxility etc., and the size of the colloids that sets the scale for variations of the elastic degrees of freedom, i.e., the nematic director field, $`𝐧(𝐫)`$ . Both microscopic (molecular dynamics simulations , tensor order-parameter free energy ) and macroscopic (elastic free energy minimization using: trial functions , simulated anealing , finite elements  and electrostatic analogies ) approaches have been used to calculate the orientational order of a nematic in the presence of a single colloid, in 2 and 3 dimensions.

Theoretical  and experimental  work on inverted nematic emulsions in 3D, reveals that for sufficiently strong homeotropic anchoring, spherical droplets induce topological defects, i.e., singularities in the nematic director field. There are several of these defects and their stability depends on parameters such as the size of the droplet, the anchoring strength and the boundary conditions. When the size of the droplet is large, a hyperbolic point defect appears at a certain distance from the droplet, and this structure is known as the satellite configuration . For smaller droplets and/or weaker anchoring strengths, the minimum elastic free energy corresponds to a ring disclination that surrounds the droplet on the equatorial plane, and this structure is known as the saturn-ring configuration . The distortion of the director field due to the presence of more than one of these colloids induces multipolar long-range interactions between the colloids . Depending on the symmetry of the distortion, the dominant long-range interactions are dipolar  or quadrupolar . At short range, however, the colloids repel each other  and this prevents their colapse.

A 2D system was also investigated by Pettey et al.  who obtained analytic solutions for the smectic-$`C`$ director around circular disks. For a disk at the center of a smectic-$`C`$ layer, the solution that minimizes the elastic free energy exhibits a hyperbolic point defect, at a distance $`r_d^0=\sqrt{2}a`$ from the center. The free energy of a system of multiple disk-defect pairs was also obtained, for pairs that are far apart. For defects aligned with the far field director, a dipole-dipole interaction results between the colloids . These predictions were tested in a recent experiment by Cluzeau et al. . Thin films of smectic-$`C^{}`$ were heated and the temperature was stabilized when cholesteric ($`N^{}`$) droplets started to nucleate. Cluzeau and collaborators observed that a hyperbolic point defect is nucleated with each cholesteric droplet, at the distance predicted by theory. In addition, they observed that the droplets interact through a long-range attraction and assemble into chains, confirming the dipolar nature of the interaction. Finally, they measured a mean separation $`D=(1.3\pm 0.1)d_a`$ between nearest-neighbor droplets in chains, where $`d_a`$ is the droplet average diameter.

Even in 2D exact solutions are not readily obtained. The analytic solutions of  are restricted to one inclusion or to multiple inclusions that are sufficiently far apart. Recently, Fukuda and Yokoyama solved the 2D problem using the tensor order parameter formalism and an efficient adaptive grid method to integrate numerically the dynamic equation . One of the advantages of this method is that no special treatment of topological defects is required, by contrast with the director description, where defects appear as singularities of the director field . For all the parameters used in , the authors found that the equilibrium director configuration around a circular disk, exhibits quadrupolar symmetry and corresponds to two defects with topological charge -1/2. This contrasts with the 3D result, where the quadrupolar defect was found to be stable over a limited range of parameters. It is also in apparent contradiction with the experiment on 2D smectic-$`C^{}`$ films . The latter, however, do not exhibit mirror symmetry about reflection of the 2D director, excluding a distortion with quadrupolar symmetry .

In this paper we use finite elements and adaptive meshes to minimize the 2D elastic free energy functional. We focus on dipolar distortion fields as observed in the experiments . In particular, we study the interaction between two disks of equal radii, $`a`$, in a large smectic-$`C`$ layer. The minimal free energy occurs when the two dipoles are aligned head to tail, in a direction parallel to the far field director. In this geometry we investigated a wide range of disk separations. We verified the validity of the electrostatic analogy at large distances and calculated the short-range repulsive interaction between two disks.

## II Model and numerical method

The free energy of a nematic is invariant under uniform rotations of the whole sample and under the symmetry operations $`𝐧𝐧`$ and $`𝐫𝐫`$. In addition, $`𝐧`$ is a unit vector so that $`n_i_jn_i`$ is zero. These considerations imply that the stiffness tensor has 2 independent bulk components in 2D. The distortions whose energy is measured by these constants are (1) splay, with non-zero $`𝐧`$, and (2) bend, with non-zero $`\times 𝐧`$. The elastic free energy for the bulk nematic phase, is the sum of these terms :

$$F=\frac{1}{2}[K_S(𝐧)^2+K_B(\times 𝐧)^2]𝑑S$$
(1)
where $`K_S`$ and $`K_B`$ are the elastic constants associated with splay and bend distortions, and $`𝐧=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$ is the 2D nematic director. The integral is over the region occupied by the 2D nematic film.

Note that if $`K_S=K_B=K`$, the elastic free energy is proportional to $`_{ij}_in_j_in_j`$ (plus boundary terms) which is the form produced by simple models of the nematic state, such as those with pairwise intermolecular interactions depending only on the angle between the long molecular axes or the 2D XY-model .

$`F`$ $`=`$ $`{\displaystyle \frac{K}{2}}{\displaystyle [(𝐧)^2+(\times 𝐧)^2]𝑑S}`$ (2)
$`=`$ $`{\displaystyle \frac{K}{2}}{\displaystyle (\theta )^2𝑑S}`$ (3)
Line and boundary terms have been ignored in (3) which is justified in the strong anchoring regime (fixed orientation at the boundaries) for systems with large line tensions (inclusions with fixed shape) since their contribution to the free energy is constant. The one elastic constant approximation (3) is not strictly necessary, since the numerical method to be described below can be applied to the Frank free energy (1); the choice of (3) however, allows a direct comparison with the analytical results of  and provides a stringent test of our numerics.

The geometry used in the numerical calculations is shown schematically in figure 1a. We consider two circular disks of radius $`a`$ separated by a distance $`D`$. Close to each disk we pin a defect at distances $`r_{d_1}`$ and $`r_{d_2}`$ measured from the center of each disk, respectively. The system has mirror symmetry with respect to $`x`$-axis. Thus, the center of the disks and the accompanying defects are located on the $`x`$-axis. The smaller circles around the defects are the defect cores (regions where the nematic order is destroyed ). They are used as cut-offs for the integral of the elastic free energy density (3). Lengths are measured in units of $`a`$ and energy is measured in units of $`K`$.

We fixed the nematic orientation, $`\theta `$, at the physical boundaries. The far field director was taken parallel to the $`x`$ axis ($`\theta =0`$) and homeotropic boundary conditions were imposed at the disk perimeters. The solution for the director field is found by minimizing the total free energy of the system, using standard numerical procedures . The major difficulty in the numerical problem stems from the large difference in the length scales set by the disk and the defect. It was overcome by using adaptive meshing techniques. With the help of a 2D mesh generator  (BL2D package), a first triangulation respecting the predefined physical boundaries, as well as the defect cores, is constructed. A final grid for the region close to a defect is shown in figure 1b and was obtained after several adaptation loops. Note that in addition to the outer boundaries (figure 1a), the triangulation follows the perimeter of the defect core. Using standard numerical procedures the elastic free energy is minimized in a region that includes the defect core. The elastic free energy of the core region is subtracted at the end of the calculation. The core energy was set to zero .

Adapted meshes may be generated if the Hessian of the solution, $`^2\theta /xy`$, is known. Far from the disks, $`\theta `$ varies slowly and the triangles can be large. By contrast, close to the defects, the Hessian is large (it diverges at $`r_d`$) and the triangles have to be made several orders of magnitude smaller. The Hessian, $`_{ij}^k`$, at vertex $`k`$ is estimated  by using the following ’weak’ definition (, p. 349)

$$_{ij}^k=\frac{\frac{\theta }{x_i}\frac{v^k}{x_j}𝑑S}{v^k𝑑S}.$$
(4)
where the integral is over the region occupied by the smectic layer. $`v^k`$ is the piecewise linear hat function associated with vertex $`k`$ ($`v^k=1`$ at vertex $`k`$ and $`v^k=0`$ elsewhere), $`x_1=x`$ and $`x_2=y`$. To construct an anisotropic non-uniform mesh, we require a metric map that distributes the interpolation error, i.e., the difference between the interpolated function and the exact solution, in a uniform fashion . The required map, $`_{ij}^k`$, at vertex $`k`$ is a two dimensional positive definite matrix given by

$$^k=c_0𝒪\left(\begin{array}{cc}|\lambda _1|& 0\\ 0& |\lambda _2|\end{array}\right)𝒪^1$$
(5)
where $`𝒪`$ is the orthogonal matrix that diagonalizes the Hessian $`^k`$, and $`\lambda _1`$ and $`\lambda _2`$ the corresponding eigenvalues. The constant $`c_0`$ controls the number of mesh triangles. To avoid problems caused by the divergence of the Hessian, the eigenvalues are bounded by

$$\lambda _i=\mathrm{min}(\mathrm{max}(|\lambda _i|,\frac{1}{c_0h_{max}^2}),\frac{1}{c_0h_{min}^2})$$
(6)
where $`h_{min}^2`$ and $`h_{max}^2`$ are the minimal and maximal mesh edge lengths.

As a test of the numerics, we considered one single disk in a smectic-$`C`$ film. We calculated the free energy as a function of the defect distance $`r_d`$, for a dipole parallel to the far field director. The equilibrium free energy was found for a defect at $`r_d=1.41\pm 0.01`$, in good agreement with the analytical result for small defects ($`r_d=\sqrt{2}`$)  and in line with the experimental result for inclusions of cholesteric droplets in smectic-$`C^{}`$ films ($`r_d=1.4\pm 0.1`$) .

## III Results

In the calculations, we used a rectangular integration region of size $`L_x\times L_y=50\times 20`$ (see figure 1a). With the adaptive mesh technique we have studied defects with cores as small as $`10^3`$, with a minimal mesh edge length of $`10^4`$. During the adaptation loops this limit was reached close to and inside the cores, providing a good discretization of the core region (see figure 1b). The final mesh configuration contained up to $`10^4`$ vertices. After 4 or 5 iterations of the adaptive mesh procedure, the relative error of the free energy is of the order of $`10^4`$ (see figure 1c). For a particular geometry, the central part of the final mesh is shown in figure 2a. In figure 2b, a gray scale was used to plot the corresponding equilibrium tilt angle, $`\theta (x,y)`$. White is $`\theta =0`$ and black $`\theta =\pi `$.

Near the core (in a small region of size $`<10^2`$) the elastic distortion field is determined by the defect alone. We have checked numerically that a change in the core size from $`\xi _1`$ to $`\xi _2`$ leads to a free energy difference of $`\pi K\mathrm{log}(\xi _1/\xi _2)`$ , independent of the position of the disks.

For these systems (small defect cores) we found that the separation of the disks had no influence on the position of the outer defect, $`r_{d_2}`$ : the minimum free energy was obtained for an outer defect at $`r_d^0`$, the position of the defect for an isolated disk, i.e., independent of the presence of the first disk. In what follows we set $`r_{d_2}=1.41`$.

The elastic free energy as a function of the position of the inner defect, $`r_{d_1}`$, is shown in figure 3a, for several disk separations. The lowest free energy was found for a disk separation $`D=2.82`$, with an inner defect at $`r_{d_1}=1.41`$.

In figure 3b we plot the position of the inner defect, $`r_{d_1}`$, that minimizes the elastic free energy, as a function of the disk separation, $`D`$. We identify three regimes. For small disk separations ($`2<D2.82`$) the position of the inner defect, $`r_{d_1}`$, is symmetrical with respect to the disks and scales linearly with $`D`$, namely $`r_{d_1}=D/2`$, i.e., the defect is at the mid-point between the disks. In the intermediate regime ($`2.82<D<6.5`$) the defect position varies non-monotonically with D, reaching a maximum, $`r_{d_1}=1.52`$, at $`D=3.5`$. Finally, for large disk separations ($`D>6.5`$) the position of the inner defect is independent of D and is given by that of an isolated disk, $`r_{d_1}=1.41`$, i.e., it is independent of the presence of the second disk.

Note that the minimal free energy corresponds to a symmetrical configuration: the inner defect is at a distance from the first disk identical to the distance of the outer defect from the second disk.

Finally, we plot the elastic free energy as a function of the disk separation, $`D`$, in figure 4a. The free energy exhibits a pronounced minimum at $`D=2.82\pm 0.01`$, close to the value $`2.6\pm 0.2`$ observed experimentally by Cluzeau et al. , for chains of (polydispersed) circular droplets. Note that the experimental estimate of D, if taken without error bars, suggests that the inner defect is nucleated at a distance smaller than the distance corresponding to isolated droplet-defect pairs. This configuration, however, is ruled out by our model (see figure 3a, open circles and squares). Within the (large) error bars we conclude that the experimental result is in line with the theoretical prediction: the lowest free energy corresponds to a symmetrical configuration, with the inner defect at a distance from the disks that is identical to the distance of deffects in isolated disk-defect pairs.

Admitedly, some effects were left out from our calculation, the most important of which is arguably that of fluctuations. In fact, fluctuations were predicted to be important in 2D  but were reported to be negligible in the rather thick films studied by Cluzeau et al. . It is possible that many-body interactions may account for a reduction in disk separation for disks in chains, when compared to the equilibrium separation of a pair of disks. This could be checked within our model by calculating the free energy of three interacting disks. The experimental results, however, suggest that this effect, if present, is of the order of the experimental error. A smaller disk separation, than that predicted by our model, may also result from surface contributions that were oversimplified in (3). In fact, deviations from strong homeotropic anchoring at the disks reduce the effective size of the inclusions, leading to a smaller equilibrium distance between disks .

The signature of the dipolar interaction for large disk separations may be seen in figure 4b. Here, the circles are the (calculated) force between disks, as a function of their separation. The force is clearly a power law that decays as $`D^3`$. Note that the numerical results for the force are less accurate than those of the free energy, since they were obtained by numerical differentiation of the latter, figure 4a.

## IV Summary and Conclusion

We studied the interaction of two circular disks immersed in smectic-$`C`$ films. For strong homeotropic anchoring the circular disks nucleate topological defects of the 2D director field. We implemented a finite elements, with adaptive meshing, numerical method capable of dealing with the variation of the director field over the widely different length scales set by the disks and the defects. The method allows the construction of well defined core regions around the defects and leads to very accurate results for the elastic free energy.

We have verified the long-range dipolar character of the interaction between disk-defect pairs, predicted by Pettey et al. . In addition, we have calculated the interaction at small disk separations, and found that at static equilibrium, the distance between the centers of the disks is equal to $`D=2.82\pm 0.01`$.

Finally, we have determined the position of the inner defect, as a function of the separation between disks. For small disk separations the defect is located symmetrically at the mid-point between the disks. At intermediate separations, this symmetry is broken and the position of the defect reaches a maximum $`r_{d_1}=1.52`$ for a disk separation, $`D=3.5`$. Finally, at large disk separations the defect relaxes slowly to the equilibrium position, $`r_d^0`$, of an isolated disk.

We intend to use the numerical method described in this paper to study the properties of other reverse nematic emulsions in 2 and 3D. In particular, we will consider the quadrupolar defects, recently studied by Fukuda and Yokoyama  using a tensor order parameter formalism. Although these authors concluded that the dipolar configuration is dynamically unstable, and the quadrupolar configuration is the only one allowed in 2D nematics, the energetics of these two types of pinned defects has not been compared.

Acknowledgments

PP thanks J. M. Tavares for get him interested in this problem and for helpful discussions. We also thank Francisco de los Santos for useful comments.

We acknowledge the support of the Fundação para a Ciência e a Tecnologia (FCT) through a running grant (Programa Plurianual) and grants No. PRAXIS XXI/BPD/16309/98 (PP) and No. SFRH/BPD/1599/2000 (M.T.).