text:
# 1 Introduction

## 1 Introduction

Both understanding of the processes occurring in a nucleus and precise calculation of important for practice parameters of these processes (for example, neutron cross-sections) requires their detailed experimental study in the wide excitation energy interval and in the large mass range of nuclei.

In these frameworks, the two-step $`\gamma `$-cascades following thermal neutron capture in target nucleus $`{}_{}{}^{117}Sn`$ have been studied. The energies of $`\gamma `$-transitions and intensities of most strong 454 cascades exciting 303 levels in the compound nucleus were determined in the energy interval practically up to $`B_n`$.

## 2 Experiment

The experiment was performed at LWR-15 reactor in Řež . The target, enriched in $`{}_{}{}^{117}Sn`$ to more than 91%, was used in the experiment. The contribution of other $`Sn`$ isotopes in the background of $`\gamma \gamma `$ coincidences was very small.

The sum coincidence spectrum of the two-step cascades measured in the experiment is presented in Fig. 1. Using this spectrum, intensity distributions for the cascades terminating at 8 low-lying levels with $`E_f<2.8`$ MeV of $`{}_{}{}^{118}Sn`$ were built. Each two-step cascade in such spectrum is presented by a pair of peaks with equal  areas and widths. Energy resolutions in each of these distributions is varying from 1.9 keV at their ends to 3 keV in centre. Such a resolution at the achieved statistics of coincidences, low and practically constant at any energy background in these spectra, which is provided the detection threshold $`0.6\times 10^4`$ events per decay for the cascades to the ground state and $`2÷5\times 10^4`$ events per decay for the cascades at the higher-lying levels.

The spectrum in Fig. 1 is very specific as compared with those for other nuclei, studied by us earlier. It has single-escape and double-escape peaks corresponding not only to one cascade transition but also to both transitions simultaneously. Thus, cascades to the ground and first excited states of $`{}_{}{}^{118}Sn`$ create in the sum coincidence spectrum four well-expressed additional peaks shifted by 511 keV.

## 3 Decay scheme

The confidence of the data on peculiarities of the nuclear excited states determining intensities of the observed two-step cascades is provided, first of all, by the reliability and completeness of the obtained scheme of $`\gamma `$-decay. Two-step cascades following neutron capture unambiguously set energies of the initial and final levels but do not contain direct information on the energy of their intermediate level. However, if cascade is intense enough ($`i_{\gamma \gamma }10^4`$ events per decay) then it is usually observed as a pair of the resolved peaks. In this case one can determine the quanta ordering in cascade.

The method to construct a decay scheme using obvious thesis about the constancy of the energy $`E_1=B_nE_i`$ of the primary transition in the cascades with the different total energy $`E_1+E_2=B_nE_f=const`$ was described for the first time in . The method uses the multi-dimensional Gauss distribution in the framework of the maximum likelihood method in order to select probable $`\gamma `$-transitions with equal energy in different spectra. As it was shown in , the algorithm gives reliable results even at the mean error $`1`$ keV in determination of the transition energies for some hundreds of cascades placed into the decay scheme. Corresponding decay scheme for $`{}_{}{}^{118}Sn`$ is given in Table 1.

Relative intensities of most intense cascades for all 8 spectra obtained in the experiment were transformed into absolute values by normalization to the absolute intensities of their primary transitions  (Table 2) and branching ratios obtained in traditional manner from the mass of coincidences accumulated in the experiment. The involving of the maximal number of the most intense cascades in normalization procedure decreases its uncertainty due to correlations between them. Some notion of the values of statistical and systematic errors is given by a comparison of the sum of intensities of all cascades proceeding via the same intermediate level with the experimental intensity of their primary transitions. At the absence of systematic errors in determination of cascade intensities, the ratio $`R=i_{\gamma \gamma }/i_1`$ should decrease when $`E_1`$ decreases. Exceeding of $`R`$ over the unity gives some total uncertainties for both intensities. Besides, owing to sufficiently different backgrounds in the spectra of cascade transitions and single $`\gamma `$-transitions, contribution of systematic error of $`i_1`$ in this uncertainty is larger than that of $`i_{\gamma \gamma }`$ at least for the smallest $`E_1`$.

The total absolute intensities $`I_{\gamma \gamma }=i_{\gamma \gamma }`$ of cascades with a fixed sum energy (including those unresolved experimentally) obtained in this way at the detection threshold of quanta $`E_\gamma >520`$ keV are listed in Table 3. Here energy interval of the cascade final levels $`E_f`$ is limited by conditions of the experiment – at the higher energies $`E_f`$ the ratio peak/background in the sum coincidence spectrum decreases. This does not allow one to get reliable information on the function $`i_{\gamma \gamma }=f(E_1)`$ for $`E_f>23`$ MeV without the Compton-suppressed spectrometer.

The total intensities $`I_{\gamma \gamma }`$ are suitable for revealing the serious discrepancies between the experiment and model notions of level density and radiative strength functions. Results of calculation of $`I_{\gamma \gamma }^{mod}`$ with the use of level density model , model of radiative strength functions  for E1- transitions and strength function k(M1)=const are listed in Table 3. Considerable exceeding of the experimental cascade intensities over the calculated values means that the level density really excited after the thermal neutron capture is less than predictions of model  and, probably, that the radiative strength functions of the primary transitions have more strong energy dependence than that predicted by model .

### 3.1 Comparison with a known decay scheme

The decay scheme of $`{}_{}{}^{118}Sn`$ includes 454 two-step cascades proceeding via 303 intermediate levels with the excitation energy up to $`(B_n1)`$ MeV; 110 levels from them are depopulated by two or more secondary transitions. Hence, reality of these levels is confirmed with a high confidence. (The mean uncertainty in determination of the energy of the cascade quantum is 0.48 keV). Quanta ordering for 191 cascades cannot be determined within algorithm  because corresponding intermediate levels are depopulated only by one secondary transition. In general case, many of these cascades can have primary transition with lower energy than energy of the secondary transition. Quanta ordering in these 191 cascades was determined using the data ENDS-file  on the estimated decay scheme of the nucleus under study.

It should be noted that the levels 2120, 2577 and 2725 keV listed in this file were not observed in our experiment neither in the sum coincidence spectrum nor as the intermediate levels of other cascades. Therefore, using data  one should take into account that the estimated decay scheme can contain wrong levels with mistaken spin values at this and, moreover, at higher energy.

Neutron capture cross section in a target is mainly determined by the only neutron resonance with the known parameters. Therefore, one can unambiguously assign spin and parity $`J^\pi =1^+`$ to the compound state of $`{}_{}{}^{118}Sn`$. So, levels of this nucleus excited by the cascades (see Table 1) have, most probably, both positive and negative parity and spin values in the interval 0-2. The absence in the sum coincidence spectrum of the peak corresponding to registration of the two-step cascades to level $`E_f=2280`$ keV, $`J^\pi =4^+`$  shows that, as in nuclei studied earlier, the total intensity of cascades with one purely quadrupole transition in $`{}_{}{}^{118}Sn`$ is less than the detection threshold of the experiment.

### 3.2 Fluctuations of the cascade intensities and problem of completeness of the decay scheme

According to the conventional notions, partial widths of primary transitions following decay of neutron resonance of non-magic nucleus are the random values whose fluctuations with respect to the average obey the Porter-Thomas distribution . Unfortunately, this notion was not verified for both primary $`\gamma `$-transitions to the levels with the energy about several MeV and in the region of the smallest partial widths. Hence, non-Gaussian character of the distribution of some amplitudes of the primary $`\gamma `$-transitions of a given multipolarity cannot be excluded. This can be due to, for example, contribution of large components of the wave function of the levels connected by corresponding transition in the matrix element of the primary $`\gamma `$-transition. In this case, experimental distribution of small partial widths will have less probability than it follows from model  owing to violation of the condition of application of the central limit theorem of the mathematical statistics.

The necessity of this question in practical solution arises at the analysis of the experimental data on the cascade intensities. Very low and rather even background in intensity distributions of the cascades to the lowest levels of the studied nucleus provides better sensitivity for the experiment on the study of coincidences relative to the traditional experiment  with a single detector. This problem for $`{}_{}{}^{118}Sn`$ is more important than, for example, for deformed nuclei due to bigger number of intermediate levels, which are depopulated by the only $`\gamma `$-transition with the intensity exceeding registration threshold.

The presence of 191 cascades of this kind and, correspondingly, the necessity to determine energies of their intermediate levels (without using method ) requires one to solve the dilemma:

(a) either there are many levels at relatively low excitation energies in $`{}_{}{}^{118}Sn`$ (for instance, $`E_i<34`$ MeV), which do not manifest themselves in other reactions, and, as a consequence, level density in this nucleus is noticeably higher than even the predictions of model , at least, in this energy interval;

(b) or the number of primary transitions whose widths are 5-10 times smaller than the mean value is considerably less than expectation according to  for the same excitation energy interval.

Considerably less threshold value of the function approximating the distribution for $`E_m>5.5`$ MeV in comparison with this parameter for lower energies admits only two these possible interpretations. This conclusion follows from Fig. 2, which represents the comparison between the cumulative sums of the experimental intensities of cascades proceeding via the same intermediate level (for the energy bin $`\mathrm{\Delta }E_i=0.5`$ MeV) and corresponding approximation (as it was done for the first time in ). Experimental data in Fig. 2 were obtained under condition that the low-intensity cascades were placed into the decay scheme in such a manner that their primary transitions have the lower energy. This assumption was used only when it was impossible to use for this aim the method  or data on excited levels from .

## 4 Cascades with the low-energy primary $`\gamma `$-transitions

Fig. 3 represents a part of the experimental intensity distribution of the cascades to the ground state of $`{}_{}{}^{118}Sn`$. Decay scheme of this nucleus up to $`E_{ex}2.5`$ MeV was very reliably established in different experiments. Based on all available data one can state to a precision of about 100% that this spectrum contains in the interval 520 to 2500 keV only 4 peaks corresponding to registration of the low-energy secondary transitions  of cascades. The rest of the spectrum should be related with the registration of the cascades with the energies of primary transitions $`E_1<2.5`$ MeV. The majority of the cascades with the high-energy primary transitions are usually registered in resolved intense peaks, and a number of low-intensity cascades with low-energy primary transitions form pseudo-continuous distribution in the spectrum.

But observation of the resolved peak related with registration of possible low-energy primary transition of corresponding cascade permits one to get some information on structure of intermediate level of cascade. Intensity of individual cascade proceeding via a given intermediate level:

$$i_{\gamma \gamma }=i_1\times i_2/i_2$$
(1)
is always less than the intensity of its primary transition $`i_1`$ because the branching ratio $`b_r=i_2/i_2`$ cannot be greater than the unity. Unfortunately, background under the peaks in the sum coincidence spectrum increases when the total energy of cascade $`E_1+E_2`$ decreases. This makes impossible the extraction of the main part of the secondary transition intensities at the decay of any high-energy enough intermediate level and, as a consequence, to determine $`i_1`$ as a sum of all values $`i_{\gamma \gamma }`$ for a given intermediate level. But for some probable primary transitions one can select from  close on the energy $`\gamma `$-quantum (which in many cases was not placed into the decay scheme suggested by authors or placed in such a manner that it was not observed in corresponding distribution of the cascade intensities). Radiative strength functions determined for such cases can be compared with the predictions of model . Analysis , made for the experimental intensity distributions $`I_{\gamma \gamma }`$ of two-step cascades to the ground and first excited states of $`{}_{}{}^{118}Sn`$ with accounting for the known total radiative width $`\mathrm{\Gamma }_\lambda =80`$ meV, shows that the radiative strength function up to the energy $`E_13.5`$ MeV does not exceed the value predicted by model . Therefore, this model can be used for estimation of probability of low-energy primary transitions whose partial widths considerably exceed the mean value.

In this experiment have been observed, at least, 15 cascades with the most intense, low-energy primary transitions ($`E_1<2.5`$ MeV) whose partial widths are, at least, 20 times larger than the prediction of . According to model , approximately 26000 levels are expected to be excited by the primary $`E1`$ and $`M1`$ $`\gamma `$-transitions in the excitation energy region from 6.5 to 8.3 MeV. The probability to observe the only random exceeding of partial width over the average by a factor of 20 and more for the considered $`\gamma `$-transitions in the framework of the Porter-Thomas distribution is equal to $`10^5`$.

Hence, from the above statement that some of the observed in  $`\gamma `$-transitions are really primary follows that some part of these transitions have anomalous large widths whose appearance cannot be random. Therefore, some of the levels with the excitation energy $`6`$ MeV must considerably differ in the structure of the wave functions from the neighbouring levels. I. e., the effects like those observed at the low excitation energy manifest themselves and at higher excitations.

The presence of very intense low-energy primary $`\gamma `$-transitions was established earlier  in the $`{}_{}{}^{126,127,129,131}Te`$ compound nuclei. They are interpreted as direct capture of neutron at the sub-shells $`3P_{3/2}`$ and $`3P_{1/2}`$. Unfortunately, clarifying of possible structure of the considered levels of $`{}_{}{}^{118}Sn`$ is impossible without further experimental and theoretical investigations.

## 5 Possible regularity of the excitation spectrum of the intermediate levels of the most intense cascades

According to the modern theoretical notions, the wave function structure of any excited states is determined by co-existence and interaction between the fermion (quasi-particles) and boson (phonons) excitations. With the excitation energy a nucleus transits from practically mono-component excitations of the mentioned types to the mixed (quasi-particles $``$ phonons) states with rather different  degree of their fragmentation. This process should be investigated in details but there is no adequate experimental methods to study the structure of the wave functions, at least, up to the excitation energy not less than 5 MeV.

Nevertheless, some information on the probable dominant components of wave functions of heavy deformed nuclei can be obtained even in this case. The authors of  suggested to search for the regularity in the excitation spectra of the intermediate levels of the most intense cascades by means of auto-correlation analysis of the smoothed distributions of the sum cascade intensities from Table 1. Intensities were smoothed by means of the Gaussian function: $`F(E)=i_{\gamma \gamma }\times exp(0.5(\mathrm{\Delta }E/\sigma )^2)`$. The distribution of this type smoothed with the parameter $`\sigma =25`$ keV is given in Fig. 4, and the values of the auto-correlation function

$$A(T)=\underset{E}{}F(E)\times F(E+T)\times F(E+2T)$$
(2)
for different selection thresholds of intense cascades are shown in Fig. 5. As it was shown in , such an analysis cannot give unique value of the equidistant period $`T`$ even for the simulated spectra (for example, for 25 “bands” consisting from 4 levels with slightly distorted equidistant period and intensities of cascades) and provide estimation of the confidence level of the observed effect. In principle, both problems can be solved in the experiments on the study of the two-step cascades in different resonances of the same nucleus. But some grounds to state that the regularity really exists can be obtained from a comparison of the most probable equidistant periods in different nuclei. The set of the probable equidistant periods obtained so far (Fig. 6) allows an assumption that the $`T`$ value is approximately proportional to the number of boson pairs of the unfilled nuclear shells. This allows one to consider the effect at the level of working hypothesis, though the probability of random existence of the regular spectrum of intermediate levels of the most intense cascades in each, arbitrary nucleus cannot be equal to zero.

The regularity in the excitation spectra testifies to the harmonic nuclear vibrations. Thus, one can assume that the structure of the intermediate levels of the studied cascades contains considerable components of the rather weakly fragmented states like multi-quasi-particle excitations $``$ phonon or several phonons. This provides logical explanation of serious decrease in the observed level density as compared with the predictions of the non-interacting Fermi-gas model: nuclear excitation energy concentrates on phonos but quasi-particle states in the energy interval up to $`2`$ MeV are excited weakly or very weakly due to insufficiency of energy for breaking of paired nucleons. This is quite probable but, in principle, not the only explanation.

## 6 Cascade population of levels in $`{}_{}{}^{118}Sn`$.

Obvious deviations of level density and radiative strength functions from predictions of from the simplest models like  and probable regularity in the excitation spectra of intermediate levels of the most intense cascades show that the notions of the cascade $`\gamma `$-decay process need development and some corrections. This necessity follows from a comparison between the experimental and model calculated values of the population of the high-lying levels of the nucleus under study. It is seen from eq. (1) that at the presence of reliable information on the intensities of individual cascades and known  intensities of their primary and secondary transitions one can to determine $`i_2`$, which equals the total population of a given level by the direct primary transitions and any possible cascades. It is worthwhile to subtract the primary transition intensity from this value in order to decrease its fluctuation. The populations obtained in this manner for the levels with positive parity are compared with the model calculated values in Fig. 7. It should be noted that difference between the calculated populations of levels with the same spin but different parities is insignificant. It is seen that the correspondence between the experiment and calculation for the levels above $`3`$ MeV cannot be achieved. On the one side, this result confirms the conclusion made above on the presence of the excitation energy interval where occurs principle change in properties of the excited states of not only deformed but also spherical nuclei and, on the other side, it points at the necessity to develop more precise nuclear models for the all interval of the excitation energies $`E_mB_n`$.

## 7 Conclusion

For the first time, the large enough and reliable scheme of the excited states and modes of their decay for the $`{}_{}{}^{118}Sn`$ compound nucleus were obtained in the experiment with the excitation energy more than 3 MeV.

Analysis of this information shows, in particular, that an improvement of precision in description of properties of spherical nucleus requires more detailed accounting for co-existence and interaction of quasi-particle and phonon components of wave functions of the levels excited after the neutron capture.

This work was supported by GACR under contract No. 202/03/0891 and by RFBR Grant No. 99-02-17863.

References

| $`[1]`$ | J. Honzátko et al., Nucl. Instr. and Meth. A376 (1996) 434 |
| --- | --- |
| $`[2]`$ | A.M. Sukhovoj and V.A. Khitrov, Sov. J.: Prib. Tekhn. Eksp. 5 (1984) 27 |
| $`[3]`$ | Yu.P. Popov, A.M. Sukhovoj, V.A. Khitrov and Yu.S. Yazvitsky, Izv. AN SSSR, Ser. Fiz. 48 (1984) 1830 |
| $`[4]`$ | S.T. Boneva, E.V. Vasilieva and A.M. Sukhovoj, Izv. RAN., Ser. Fiz., 51(11) (1989) 2023 |
| $`[5]`$ | Yu.E. Loginov, L.M. Smotritskij, P.A. Sushkov, Vopr. Atom. Nauki i Tech., ser. Fiz. Yad. Reac. (in Russian), 1/2 (2001) 72 |
| $`[6]`$ | W. Dilg, W. Schantl, H. Vonach and M. Uhl, Nucl. Phys. A217 (1973) 269 |
| $`[7]`$ | S.G. Kadmenskij, V.P. Markushev and W.I. Furman, Sov. J. Nucl. Phys. 37 (1983) 165 |
| $`[8]`$ | K. Kitao, Nucl. Data Sheets, 75(1) (1995) 99 |
|  | http://www.nndc.bnl.gov/nndc/nsr/nsrframe.html |
| $`[9]`$ | C.F. Porter and R.G. Thomas, Phys. Rev. 1956. 104, (1956) 483 |
| $`[10]`$ | A.M.Sukhovoj and V.A. Khitrov, Physics of Atomic Nuclei, 62(1) (1999) 19 |
| $`[11]`$ | E.V. Vasilieva, A.M.Sukhovoj and V.A. Khitrov, Physics of Atomic Nuclei, 64(2) (2001) 153 |
| $`[12]`$ | I.Tomandl, V.Bondarenko, D.Bucurescu, J.Honzátko, T.von Egidy, H.-F.Wirth, G.Graw, R. Hertenberger, A.Metz, Y.Eiserman, Proc. 10 International Symposium on Capture gamma-ray spectroscopy and related topics, Santa Fe, USA, Ed. by S.Wender, AIP, (1999) 200 |
| $`[13]`$ | L.A. Malov and V.G. Soloviev, Yad. Fiz., 26(4) (1977) 729 |
| $`[14]`$ | A.M. Sukhovoj and V.A. Khitrov, Bulletin of the Russian Academy of Sciences, Physics 61(11) (1997) 1611 |
| $`[15]`$ | E.V. Vasilieva et al., Bulletin of the Russian Academy of Sciences, Physics 57 (1993) 1582 |

Table 2.
The energies $`E_\gamma `$ and absolute intensities (% per decay) $`i_1`$ of the most intense transitions used for normalization of the cascade intensities in $`{}_{}{}^{118}Sn`$. $`i_{\gamma \gamma }`$ is the observed intensity of the cascades with corresponding primary transition.

| $`E_1`$, keV | $`i_1`$  | $`i_{\gamma \gamma }`$ |
| --- | --- | --- |
| 8096.37(17) | 0.23(3) | 0.37(2) |
| 7268.84(13) | 0.40(3) | 0.44(1) |
| 6997.7(4) | 0.37(8) | 0.31(2) |
| 6648.52(7) | 1.00(6) | 0.75(3) |
| 6422.14(11) | 1.04(9) | 0.88(4) |
| 6110.0(2) | 0.99(2) | 0.80(3) |
| 6055.66(18) | 0.52(2) | 0.55(3) |
| 5444.9(3) | 0.55(7) | 0.35(1) |
| 5298.1(5) | 0.33(5) | 0.28(1) |
| 5217.51(10) | 0.95(6) | 0.60(3) |
| 5208.7(2) | 1.1(1) | 1.09(4) |
| 4781.55(16) | 0.79(6) | 0.28(2) |
| 4765.99(15) | 0.80(5) | 0.51(3) |
| 4556.31(14) | 1.33(6) | 0.39(4) |
| Sum | 9.85(21) | 7.6(1) |

Table 3. Energies $`E_1+E_2`$ of cascades and their absolute intensities $`I_{\gamma \gamma }`$ (% per decay). $`E_f`$ is the energy of the cascade final level

| $`E_1+E_2`$, keV | $`E_f`$, keV | $`I_{\gamma \gamma }^{exp}`$ | $`I_{\gamma \gamma }^{mod}`$ |
| --- | --- | --- | --- |
| 9326.30 | 0 | 16.0(34) | 6.7 |
| 8096.63 | 1230 | 15.3(11) | 7.2 |
| 7568.00 | 1758 | 2.4(7) | 1.1 |
| 7283.42 | 2042 | 3.3(16) | 2.4 |
| 7269.39 | 2057 | 2.8(9) | 0.8 |
| 7000.36 | 2325+2328 | 5.6(9) | 2.8 |
| 6923.08 | 2403 | 2.8(2) | 1.5 |
| 6829.42 | 2497 |  | 0.5 |
| 6648.95 | 2677 | \[1.5\] | 1.0 |
| 6588.29 | 2738 |  | 1.7 |
| sum |  | 55.8(43) | 25.7 |

Note: intensity of the cascades to the levels $`E_f=`$ 2497, 2677 and 2738 keV was estimated from a comparison of their peak areas in the sum coincidence spectrum with the peak area corresponding to the cascades to the levels 2325 and 2328 keV.


format: text

text:
# Theorem.

On a sequence related to the Josephus problem
Ralf Stephan***The author can be reached at mailto:ralf@ark.in-berlin.de.

In this short note, we show that an integer sequence defined on the minimum of differences between divisor complements of its partial products is connected with the Josephus problem (q=3).

We prove the following theorem and, finally, state the relatedness of two constants.

###### Theorem.

Let $`a_n`$ and $`b_n`$ be recursively defined as

$`a_0`$ $`=4,a_n=\mathrm{min}\left(|d_jp_n/d_j|>1\right),p_n=_{k=0}^{n1}a_k,`$
$`d_j|p_n,1j\sigma (p_n).`$
$`b_1`$ $`=1,b_n=\frac{1}{2}_{k=1}^{n1}b_k.`$
(1) $`\text{Then }a_n`$ $`=2^{b_n},\text{for }n>2.`$
The first terms of $`a_n`$ and $`b_n`$ are \[S\]\[Z\]

$`a_{n0}`$ $`=\{4,3,4,2,4,8,16,64,\mathrm{}\},b_{n1}=\{1,1,1,2,3,4,6,9,14\mathrm{}\}.`$
We need two lemmata.

###### Lemma 1.

For $`k>0`$,

(2)
$$\sigma (32^k)=2k+2.$$

###### Proof.

This is true for $`k=1`$, and the set of divisors of $`32^{k+1}`$ is the set of divisors of $`32^k`$ plus $`2^{k+1}`$ and $`32^{k+1}`$ itself. ∎

###### Lemma 2.

Let $`\delta (m)`$ denote the smallest absolute value of the differences between complementary divisors of $`m>1`$:

$`\delta (m)`$ $`=\mathrm{min}\left(\left|d_j{\displaystyle \frac{m}{d_j}}\right|\right),d_j|m,1j\sigma (m).`$
Then
(3) $`\delta (32^k)`$ $`=2^{k/2},k>0.`$
###### Proof.

Let us sort the divisors of $`32^k`$ by size and call these $`D_j`$:

$$D_{1j2k+2}=\{1,(2,3),\mathrm{},(2^i,\frac{3}{2}2^i),(2^{i+1},\frac{3}{2}2^{i+1}),\mathrm{},(2^k,\frac{3}{2}2^k),32^k\}.$$

Any smallest complementary divisor difference must be the one where the divisors are in the exact middle of the sorted list, which, using (2), is $`k+1`$. And so, $`\delta (32^k)=D_{k+2}D_{k+1}`$.

Now, the proposition (3) is true for $`k=1`$. For every increase of $`k`$ by one, $`\sigma (m)`$ increases by two, and the index of the wanted pair of divisors increases by one, so $`D_{k+2}D_{k+1}`$ goes through the values

$`\frac{3}{2}2^i2^i`$ $`=2^{i1}`$
$`2^{i+1}\frac{3}{2}2^i`$ $`=2^{i1}`$
$`\frac{3}{2}2^{i+1}2^{i+1}`$ $`=2^i`$
$`2^{i+2}\frac{3}{2}2^{i+1}`$ $`=2^i`$
so it doubles every second step which is just the meaning of (3). ∎

Fixing the induction base at $`\delta (p_3=48)=2^1=2^{b_3}`$ to make sure that $`D_{k+2}D_{k+1}>1`$, the main proposition (1) is now obvious, since the powers of two in $`a_n`$ behave the same way under multiplication as unity does in $`b_n`$ under addition.

Because the asymptotics of $`b_n`$ are known\[C\], with

$$b_n=c\left(\frac{3}{2}\right)^n\frac{1}{2},c=0.36050455619661495910154466\mathrm{},$$

the investigation of $`a_{n3}=2^{b_n}`$ is settled, except for the closed form for $`c`$. Reble already proved\[R\] that $`b_n`$ is connected to the Josephus problem. Independently, our numerics show that

(4)
$$c=\frac{2}{9}K(3),$$

with $`K(3)`$ the universal constant in the same problem with $`q=3`$, a constant already discussed (\[OW\]\[HH\]), and whose closed form is still unknown.


text:
# Non-spherical evolution of the line-driven wind instability

## 1 Introduction

The winds of hot, luminous, early-type stars are composed of material photoionized by the bright ultraviolet continuum of the star. Lucy & Solomon (1970) proposed that resonance line scattering by metal ions would provide a strong source of driving, to explain the high mass-loss rates and flow velocities observed. While a single resonance line of gas at rest would only be driven by a small part of the stellar spectrum, the many ($`10^3`$) resonance lines and the variation of the line frequency which results from the Doppler shift multiply the line force until it dominates other driving mechanisms.

Detailed calculations can be made in which an accurate line-list is used and radiative transfer is calculated using a photospheric spectral energy density distribution, but the computational expense required necessitates other approximations and the accuracy of the results depends on the accuracy of the line-list used. Castor, Abbott and Klein \[Castor et al. 1975\] (hereafter CAK) proposed an alternative approach, where instead of using a detailed table of line opacities, the force is calculated assuming a distribution function of the numbers of lines present as a function of the line opacity, generally specified as a power-law (indeed, the results of more detailed studies are often characterised in terms of an effective power-law parameter). CAK use the large-velocity gradient approximation (Sobolev 1960) to allow the steady-state structures of the winds to be determined by integration of a simple ordinary differential equation.

The Sobolev approximation is more suitable in a steady state model than for a time-dependent simulation, since effects that evolve on a small scale will not be adequently accommodated. This issue was addressed by Owocki, Castor and Rybicki \[Owocki et al. 1988\] (hereafter OCR) where the Sobolev approximation is replaced by an explicit calculation of the optical depth. In this paper, they use what they term a pure-absorption approximation, in which photons which interact with the gas are removed from the stellar spectrum. While this initially seems a strange approach to modelling a wind driven predominantly by resonance scattering, the frequency redistribution (in the rest frame of the star) achieved in the scattering process means that it is a reasonable first approximation, and indeed the accuracy of the results obtained have been confirmed by subsequent, more detailed, studies by the same authors.

In the present paper, we apply the pure-absorption approximation of OCR to study the local structure and stability of line-driven winds. This work complements previous studies, which have generally either used the non-Sobolev treatment pioneered by OCR \[Owocki 1991, Owocki & Puls 1996, and further developed by\], but been restricted to spherical symmetry, or sacrificed the non-local treatment of the radiation field to allow the flow to be studied in higher dimensionalities. In common limiting cases, our results agree with these other models to reasonable accuracy. When a non-local treatment of the radiation transfer is combined with integration in a two-dimensional domain, features result which are entirely novel. An initial effort to move into the two dimensional domain was considered by Owocki (1999), where they consider a multiple ray approach to the calculation of the line driving.

The plan of this paper is as follows. In Section 2 we discuss the physical basis for a line driven from the Sobolev, CAK and OCR formalisms. In Section 3 we discuss an approximate method for the solution of the line driving force, the initial conditions used and introduce three statistical parameters for use in analysing the results. In Section 4 we introduce a 1D model, which can be compared to previous calculations. The evolution of different sizes of perturbation are considered, the data is statistically analysed and the absorbtion spectrum is calculated. In Section 5 the 1D model is extended to a 2D model and the resulting structure is discussed. A summary is presented in Section 6.

## 2 Radiative force and Hydrodynamics

The absorption and scattering of radiation is the dominant process which drives hot star winds. The interaction of stellar photons with wind particles accelerate the particles away from the star. Treated in full generality, the combined problem of hydrodynamics and radiation transfer is a highly coupled integrodifferential system. However, OCR demonstrated that a pure absorption model gave a reasonable account of the local dynamics of wind instabilities in spherical symmetry. In this paper, we study the flow of gas in less symmetrical situations, and adopt the pure absorption model as a tractable first approximation.

The general expression for monochromatic optical depth in a wind is

$$\tau _\nu (r)=_0^r\rho (r^{})\kappa _\nu (r^{})𝑑r^{},$$
(1)
where $`\rho (r)`$ is the density of the gas at a distance $`r`$ from the photosphere and $`\kappa _\nu `$ is the opacity at frequency $`\nu `$, per unit mass of the wind material. For a single absorption line, this can be expressed in terms of a profile function,

$$\tau _\nu (r)=_0^r\rho \kappa _\mathrm{l}\varphi \left(x(\nu )\frac{v(r^{})}{v_{\mathrm{th}}}\right)𝑑r^{},$$
(2)
where $`x`$ is the frequency displacement from line centre for gas at rest with respect to the star, in units of the line width, $`\mathrm{\Delta }\nu `$,

$$x=\frac{\nu \nu _0}{\mathrm{\Delta }\nu },$$
(3)
and $`\mathrm{\Delta }\nu =(v_{\mathrm{th}}/c)\nu _0`$ when the line width is dominated by thermal broadening. The lines are Doppler shifted by the local velocity, $`v(r)`$ of the stellar wind.

The profile function, $`\varphi `$, gives the relative likelihood that a photon will be scattered as a function of the offset of the frequency of the photon from line-centre in the local rest frame of the gas. It is given as a function of $`x`$, normalised so that

$$_{\mathrm{}}^{\mathrm{}}\varphi (x)𝑑x=1$$
(4)
so $`\kappa _\mathrm{l}`$ is the line opacity integrated over the line profile divided by the the line width, $`\mathrm{\Delta }\nu `$ (and so for a specific state of the wind material, $`\kappa _\mathrm{l}v_{\mathrm{th}}`$ is independent of $`v_{\mathrm{th}}`$). The line profile is determined from a combination of intrinsic (quantum) line width, thermal motion and pressure broadening. For conditions which apply in a stellar wind, it may be taken as Gaussian with a thermal width, although the Lorentzian wings may be an important factor at the base of the wind \[Poe et al. 1990\]. The variation of line profile between ions is small enough to be neglected in the present work. The role of photospheric Lorentz wings in wind acceleration in B stars was investigated by Babel (1996).

Following OCR, the acceleration due to a line at a particular frequency, $`\nu _0`$, and opacity $`\kappa _\mathrm{l}`$, can be expressed as,

$$g_{\nu _0,\kappa _\mathrm{l}}=\frac{\kappa _\mathrm{l}v_{\mathrm{th}}\nu _0F_{\nu _0}}{c^2}_{\mathrm{}}^{\mathrm{}}𝑑xl_0(x)e^{\tau (x,r)}\varphi \left(x\frac{v(r)}{v_{\mathrm{th}}}\right),$$
(5)
where the factor preceding the integral,

$$g_{\mathrm{thin}}=\frac{\kappa _\mathrm{l}v_{\mathrm{th}}\nu _0F_{\nu _0}}{c^2},$$
(6)
is the line acceleration in the optically-thin limit. The function $`l_0(x)`$ is the photospheric profile of the line, so that

$$F_{\nu _0}(r)=F_{\nu _0}l_0(x)e^{\tau (x,r)}$$
(7)
gives the strength of the local radiation field at position $`r`$ and offset $`x\mathrm{\Delta }\nu `$ from the line centre $`\nu _0`$.

In a stellar wind, there are many ionic transitions with different frequencies and line strengths. Rather than treat this ensemble of lines in detail, CAK developed a statistical approach based on a power law distribution line number density as a function of opacity. OCR extended this model, by assuming that the line opacities were distributed as

$$N(\kappa )=\frac{1}{\kappa _0}\left(\frac{\kappa }{\kappa _0}\right)^{\alpha 2}\mathrm{exp}(\kappa /\kappa _{\mathrm{max}}),$$
(8)
including an exponential cutoff at large opacity. Assuming that the underlying continuum radiation is approximately constant over the line absorption band, the radiation acceleration from the ensemble of lines is given by

$`g_{\mathrm{rad}}(r)`$ $`=`$ $`{\displaystyle \frac{N_0F\mathrm{\Gamma }(\alpha )}{c}}\left({\displaystyle \frac{v_{\mathrm{th}}}{c}}\right)^\alpha `$ (9)
$`\times `$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \frac{\varphi (xv(r)/v_{\mathrm{th}})}{(\eta (x,r)+1/\kappa _{\mathrm{max}}+\varphi (x)/\sigma _c)^\alpha }},`$
where $`\eta (x,r)`$ is the ratio of the optical depth to the opacity and can be thought of as a profile weighted column density from the photosphere to a point $`r`$ in the wind,

$$\eta (x,r)=\frac{\tau (x,r)}{\kappa _\mathrm{l}}=_r_{}^r𝑑r^{}\rho (r^{})\varphi \left(x\frac{v(r^{})}{v_{\mathrm{th}}}\right).$$
(10)
As indicated by OCR, a useful way of determining the radiation force, equation (9), is to use differences of its integral, the radiation pressure

$`\rho g_{\mathrm{rad}}`$ $`=`$ $`{\displaystyle \frac{d(P_{rad})}{dz}}`$ (11)
$`P_{\mathrm{rad}}`$ $`=`$ $`{\displaystyle \frac{N_0F\mathrm{\Gamma }(\alpha )}{(1\alpha )c}}\left({\displaystyle \frac{v_{\mathrm{th}}}{c}}\right)^\alpha `$ (12)
$`\times `$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x\left[\eta (x,r)+{\displaystyle \frac{1}{\kappa _{\mathrm{max}}}}+{\displaystyle \frac{\varphi (x)}{\sigma _c}}\right]^{1\alpha }.`$
The above neglects the effects of re-emitted or scattered radiation. One important such effect is the ‘line-drag’ phenomenon, described by Lucy (1984). This can prevent small-scale velocity perturbations from growing as a result of the anisotropic scattered radiation field. While we neglect this effect in the present work, we discuss how our methods may be extended to take account of its influence.

### 2.1 The Sobolev approximations

In the limit of narrow lines, the profile function reduces to a delta function and the integral in equation (10) can found analytically, giving

$$\tau _\nu =\underset{r_i}{}\kappa _\mathrm{l}v_{\mathrm{th}}\rho (r)\left|\frac{dv}{dr}\right|^1,$$
(13)
where the sum is over the positions, $`r_i`$, at which the line of interest is Doppler shifted to frequency $`\nu `$. This will only occur once in a monotonically accelerating wind with a single line, and is often assumed to only occur once in simple treatments of more general flows. This is known as the Sobolev, or large velocity gradient, approximation.

The local radiation driving for a single optically thick line is then given by

$$g_{\mathrm{Sob}}=\left(\frac{\nu F_\nu }{cv_{\mathrm{th}}\rho }\right)\left|\frac{dv}{dr}\right|.$$
(14)
For an ensemble of lines obeying the distribution function (8), OCR find

$$g_{\mathrm{Sob}}=\frac{N_0F\mathrm{\Gamma }(\alpha )}{(1\alpha )c}\left(\frac{1}{\rho c}\left|\frac{dv}{dr}\right|\right)^\alpha \left[\frac{(1+\tau _{\mathrm{max}})^{1\alpha }1}{\tau _{\mathrm{max}}^{1\alpha }}\right],$$
(15)
where

$$\tau _{\mathrm{max}}\frac{\rho \kappa _{\mathrm{max}}v_{\mathrm{th}}}{|dv/dr|},$$
(16)
which gives the CAK acceleration $`g_{\mathrm{CAK}}`$ in the limit $`\kappa _{\mathrm{max}}\mathrm{}`$.

The ratio of the thermal velocity to the velocity gradient, $`\mathrm{}_{\mathrm{Sob}}=v_{\mathrm{th}}/|dv/dr|`$, is called the Sobolev length. This is the length scale over which the flow velocity increases by a thermal velocity. The Sobolev approximation is invalid when velocity gradients are small and cannot be used to study the dynamics of processes with length scales smaller than the Sobolev length.

## 3 Numerical method

We performed numerical simulations of line-driven winds using the code vh-1 (Blondin 1990). This code uses the Piecewise Parabolic Method (PPM) of Colella & Woodward \[Collela & Woodward 1984\] to solve the equations of hydrodynamic grid, using Lagrangian coordinates to advance the solution, which is then re-mapped onto an Eulerian grid.

We added an additional force corresponding to line driving by a radiation field. The derivation of this force term is described in detail in the following section. The gas is assumed to be isothermal. We treat the flow relatively far from the photosphere, which allows us to use plane-parallel symmetry and to ignore finite disc effects and spherical divergence, but requires that the wind velocity is everywhere supersonic.

The scaling of the problem is kept arbitrary in the present paper, since we are interested in the general behaviour of line-driven winds. As a result, our results are not restricted to particular stellar parameters. The dimensionless parameters which must be specified include the ratio of the thermal velocity to isothermal sound speed, which we choose to be $`v_{\mathrm{th}}/c_\mathrm{s}=0.5`$ \[Poe et al. 1990\]. This is the maximum value which produces a stable flow and is still physically plausible. Lower values of this ratio are inclined to create a less steep and more unstable flow.

We take the boundary conditions to be constant velocity and density at the inner boundary and free outflow at the outer boundary throughout, with reflective boundary conditions at the sides of the two-dimensional simulations. The size of the numerical grid was kept constant, $`800`$ cells for one-dimensional simulations and $`800\times 200`$ for two dimensions.

The value of the CAK power-law line-ensemble index is fixed at a typical value of $`\alpha =0.7`$ (OCR).

### 3.1 Optical depth calculation

The calculation of the radiation force, equation (12), for a particular cell reduces to differencing the values of the radiation pressure at the interfaces of the cell,

$$P_{rad}^i_{x_2}^{x_1}𝑑x\left[\eta (x,r)+\frac{1}{\kappa _{\mathrm{max}}}+\frac{\varphi (x)}{\sigma _c}\right]^{1\alpha },$$
(17)
where $`\eta (x,r)`$ is the profile weighted column density which is the amount of radiation absorbed by the wind up to the point $`r`$. The direct calculation of the double integral implicit in this equation can dominate the time for a computation.

However, since the radiative driving force is dependent not on the radiation pressure but rather on its gradient, a local approach can be taken. At any position the transmitted spectrum will only change over a small range of offsets $`x`$ around the value for the local mean flow, and the radiation driving will be obtained from the momentum of just these absorbed photons. The force calculation can therefore be speeded up considerably by storing the local scaled optical depth $`\eta (x)`$, and advancing it by a short-characteristics method. Only the value of the previous spectrum step, $`\eta _{i1}`$, and the increase across the current cell $`\mathrm{\Delta }\eta _i`$ are necessary to calculate the driving force in cell $`i`$,

$$\frac{dP_{\mathrm{rad}}}{dz}=\underset{n=x_1}{\overset{x_2}{}}[(\eta _n+1/\kappa _{\mathrm{max}})^{1\alpha }(\eta _{n1}+1/\kappa _{\mathrm{max}})^{1\alpha }]\mathrm{\Delta }x.$$
(18)
Using the quantity

$$\eta (x,r)^{}=\eta (x,r)+\frac{1}{\kappa _{\mathrm{max}}}+\frac{\varphi (x)}{\sigma _c}$$
(19)
in our numerical calculations allows the photospheric profile and opacity limit to be absorbed in the initial conditions. In fact, we also neglect the reversing layer term, since the physical domain is assumed to be far away from the photosphere.

### 3.2 Choice of profile function

In the present paper, we approximate the profile function by a top hat function. This allows the calculation of the increase in column density to further simplified, by using the analytic formulae for the convolution of this function with the column density profile within the grid cell,

$$h_i(v)=\{\begin{array}{cc}\rho _i\mathrm{\Delta }x_i/(v_\mathrm{u}v_\mathrm{d})\hfill & v_\mathrm{d}<v<v_\mathrm{u}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}$$
(20)
$$\mathrm{\Delta }\eta _i=\varphi _ih_i,$$
(21)
where we have taken the density to be constant within grid cell and the velocity vary linearly between limits $`v_\mathrm{d}`$ and $`v_\mathrm{u}`$ at its edges. While these distributions differ from the profiles assumed in the hydrodynamical solution, we are free to assume this in the spirit of operator splitting. As we shall see, the results are sensitive to the way in which the velocity limits for the cells are determined from the cell-centred values in the hydrodynamical solution.

In reality, it would be more appropriate to take the line profile function to be a thermal Doppler function or Voigt profile. However, our approximation should not have too dramatic an effect on the results. Owocki & Rybicki (1984) compare the effect of different forms of the profile function with analytical asymptotic solutions, for the ratio of the perturbed line force, $`\delta g`$ to the perturbed velocity, $`\delta v`$. They operate the line driving on a sinusoidal velocity perturbation of the form

$$\delta v(z)=\delta v_0e^{ikz},$$
(22)
where $`\delta v_0`$ is the velocity perturbation that propagates along with the mean flow of the wind, $`k`$ is the wavenumber of the perturbation, and $`z`$ is the Cartesian distance away from the photosphere. A useful quantity to determine is the ratio of the perturbation of the acceleration to the apploed velocity perturbation, $`\delta g/\delta v`$. Owocki & Rybicki (1984) show that if this quantity has a positive real part, the flow is unstable.

When the optical depth within a Sobolev length is large, $`\tau _0=\rho _0\kappa _\mathrm{l}\mathrm{}_{\mathrm{Sob}}1`$, a top hat profile gives an analytic solution matching the asymptotic forms

$$\frac{\delta g}{\delta v}\omega _0\frac{ik}{\chi _0+ik},$$
(23)
where $`\chi _0=\rho _0\kappa `$ is the line strength in the mean flow, and $`\omega _0=g_{\mathrm{thin}}/v_{\mathrm{th}}`$ is the growth rate of perturbations in the short-wavelength limit, and the growth rate decreases to zero for long-wavelength perturbations.

In the optically thin limit, $`\tau _01`$, this ratio becomes,

$$\frac{\delta g}{\delta v}\tau _0$$
(24)
which is independent of the wavenumber and solely dependent on the optical depth. This tends to zero making the perturbed acceleration negligible.

By this perturbation analysis it is clear that while the flow feature remains optically thick the flow will be stable. This is shown in equation (23) where only the phase of the initial perturbation is affected by a change in the acceleration. In the limit of short wavenumber perturbations, $`\lambda \mathrm{}_{\mathrm{Sob}}`$ the perturbed force reduces to the Sobolev approximation $`\delta g\delta v^{}`$. The evolution of an optically thin perturbation, equation (24), is unstable since a change in the velocity field produces a change in the acceleration. Since there is no wavenumber dependence of the force in equation (24) the acceleration will continue until the perturbation becomes optically thick.

The form of the profile function will have no effect on the radiation driving force on scales larger than the Sobolev length. As we find (and as has been seen in previous work), the growth of instabilities soon leads to smooth flow regions which are separated by strong shocks and high-gradient rarefactions in which the Sobolev length will not be resolved by the numerical grid. In these latter regions, discretisation errors in the hydrodynamical scheme will probably be at least as important as errors due to our approximation to the line profile function, and so it seems reasonable to use this computationally convenient form in this first treatment.

While we have argued that the top hat-profile function is adequate for use in this initial calculation, we will in future work extend this to more realistic forms. The generalisation may be achieved by representing the profile using a sum of top hat functions, or with a higher order approximation.

As a final word of justification to the use of a top-hat function in this approximate wind solution we present a comparison of our top-hat model with preliminary results using a Gaussian profile. We defer a full comparison of the numerical effects of different profile shapes to future work, so as not to detract from the main thrust of the current paper. We compare results for the different profile functions in two limits, (i) absorption due to a wind with constant velocity gradient and constant density, (ii) absorption due to a data set, representative of a wind solution containing small scale structure. As can be seen from Fig. 1 there is little discernible difference between the two approaches. These preliminary tests indicate that the use of a top-hat function is not a serious over simplification in the limits of this numerical implementation. We do acknowledge that the Gaussian form of the profile function is a large improvement over the top-hat and for future models where smaller length and velocity scales are resolved, it should be used.

### 3.3 Initial conditions

We at first used a solution of the steady-state CAK equations to provide initial conditions for our simulations. However, it was found that the change to a discrete grid resulted in significant disturances in the subsequent flow solutions. Instead, we ran the code without perturbations and took the output data set as the starting point for subsequent simulations.

We use constant velocity ($`v=1`$) and density ($`\rho =1`$) as the inner boundary condition. The radiation incident at the inside of the grid was assumed to be a unabsorbed pure continuum. Since the flow considered here is supersonic at the outset, it is not consistent to use a photospheric profile, as assumed by OCR. As a result, there is a narrow adjustment layer where the radiation field is first absorbed, but this does not change the overall mass flux through the simulation. It will, of course, be important to include the photospheric profile in future work, when we include the geometric divergence of the flow and calculate the development of the wind from the photosphere outwards.

The values of the code units used throughout this work are shown in Table 1. Although arbitrary code units are used in the numerical simulations, using typical physical values \[Owocki et al. 1988\] of parameters (listed in Table 1), the length units used here are equivalent to stellar radii.

### 3.4 Statistical descriptors of the flow

As suggested by Runacres & Owocki \[Runacres & Owocki 2002\], the time-dependent flow of the wind can usefully be characterized by three statistical parameters: the clumping factor, velocity dispersion and correlation function. These give an indication of the mean level of flow fluctuations at different points, and the form that these take. In particular, the clumping factor will be unity and the velocity dispersion will be zero in regions where the density varies little with time, and significantly larger if the density has strong transient features. The correlation coefficient indicates if the velocity and density are in (out) of phase, corresponding to a dominance of forward (reverse) shocks in the flow, strong (anti) correlation is shown when C $`+()1`$.

The clumping factor, velocity dispersion and correlation coefficient are calculated using the relations

$`f_{\mathrm{cl},i}`$ $`=`$ $`{\displaystyle \frac{\rho _i^2}{\rho _i^2}},`$ (25)
$`v_{\mathrm{disp},i}`$ $`=`$ $`\sqrt{v_i^2v_i^2},`$ (26)
$`C_{v\mathrm{log}\rho ,i}`$ $`=`$ $`{\displaystyle \frac{v_i\mathrm{log}\rho _iv_i\mathrm{log}\rho _i}{v_{\mathrm{disp},\mathrm{i}}\sqrt{(\mathrm{log}\rho _i)^2\mathrm{log}\rho _i^2}}},`$ (27)
where the angle brackets refer to a time average of flow quantities in cell $`i`$; the parameters are functions of the position at which the time average is calculated. (The values were calculated only using discrete samples in the evolution of the simulations, so in practice we also extended the average over finite spatial regions to remove spurious numerical noise.)

## 4 One-dimensional simulations

In this section, we describe one dimensional simulations using a variety of forms for the radiation driving force. We present these in order to allow comparison of our results with those of previous authors. The simulations were performed both using a local driving force, based on the CAK formalism, and using the pure absorption approximation.

Previous studies, for example by Gayley et al. \[Gayley et al. 2001\] and Proga et al. \[Proga et al. 1998\], have applied local radiation driving terms to flows in binary stellar systems and from the surface of accretion discs. We find here that the results so obtained are crucially dependent on the manner in which the flow velocity gradient is determined.

We then present the one-dimensional results of our treatment of the pure absorption approximation. The behaviour of an unperturbed flow has been studied by many different authors (e.g. OCR, Runacres & Owocki \[Runacres & Owocki 2002\]). The basic stuctures of line driven flows are rarefactions, steep reverse shocks as fast upstream material runs into slow downstream material and de-shadowing of downstream gas, where gas recieves an increase in acceleration from the radiation field as it moves with sufficient velocity so that the incident radiation no longer lies in the range that is absorbed by the upstream patch. We find good agreement with these results. We then study the development of flows with small and large perturbations.

### 4.1 CAK force law

The radiation force may be calculated as a function of purely local conditions using the CAK formalism. As described above, this leads to an acceleration term given by equation (15) in the limit $`\tau _{\mathrm{max}}\mathrm{}`$, that is

$$g_{\mathrm{CAK}}=\frac{N_0F\mathrm{\Gamma }(\alpha )}{(1\alpha )c}\left(\frac{1}{c\rho }\left|\frac{dv}{dr}\right|\right)^\alpha .$$
(28)
This has a milder dependence on velocity gradient than that given by equation (14), as a result of allowing for the effect of locally optically-thin lines.

In order to implement the CAK driving force, a form of differencing has to be chosen for the velocity gradient. We investigate here two extreme cases. The upstream form is defined as <sup>3</sup><sup>3</sup>3This uses the convention that the inner boundary is at $`n=0`$ and the outer boundary is at $`n=n_{\mathrm{max}}`$

$$\frac{\mathrm{\Delta }v_\mathrm{u}}{\mathrm{\Delta }r}\frac{v_nv_{n1}}{r_nr_{n1}},$$
(29)
and downstream form with

$$\frac{\mathrm{\Delta }v_\mathrm{d}}{\mathrm{\Delta }r}\frac{v_{n+1}v_n}{r_{n+1}r_n}.$$
(30)
The choice of differencing direction determines the direction in which information is propagated by the radiation force term. The upstream calculation implies that the force depends on the gas properties in the current and the previous cells. In the downstream case, by contrast, gas can ‘see’ material in its shadow, which does not correspond to the underlying physics of the radiation flow (at least in the absence of scattering).

In Fig. 2, we compare the results of simulations using these two forms. We find that the downstream method, equation (30), method is stable, and recovers the solution of the steady-state CAK equations with good accuracy. When the flow is highly supersonic, the steady-state momentum equation becomes

$$v\frac{dv}{dr}=g_{\mathrm{CAK}},$$
(31)
and since the mass flux $`\mathrm{\Phi }=\rho v`$ is constant, $`g_{\mathrm{CAK}}(v|dv/dr|)^\alpha `$. As a result, the velocity increases as $`v(rr_0)^{1/2}`$, which is a good approximation to the form of the solid curve in Fig. 2.

However, the flow calculated using the upstream method, equation (29), is unstable. If a stable steady flow is obtained using the downstream method it will become unstable at every point, when the method is changed to upstream calculation. This instability is not a short lived phenomenon but persists as long as the upstream method is used. This suggests that the instability is an intrinsic property of the form of radiation driving law used, rather than simple the result of initial flow perturbations and boundary effects sweeping across the grid.

This behaviour may be understood in terms of the propagation of ‘Abbott waves’. Abbott (1980) treated the propagation of flow perturbations in the Sobolev limit. In one dimension, he found 2 radiative-acoustic wave modes moving at speeds

$`v_{}`$ $`=`$ $`v{\displaystyle \frac{1}{2}}f_\mathrm{L}^{}\sqrt{\left({\displaystyle \frac{1}{2}}f_\mathrm{L}^{}\right)^2+c_\mathrm{s}^2}`$ (32)
$`v_+`$ $`=`$ $`v{\displaystyle \frac{1}{2}}f_\mathrm{L}^{}+\sqrt{\left({\displaystyle \frac{1}{2}}f_\mathrm{L}^{}\right)^2+c_\mathrm{s}^2},`$ (33)
where $`f_\mathrm{L}^{}`$ is the derivative of the radiation force with respect to the velocity gradient. The critical points in the steady-state CAK equation occur where one or other of these mode velocities is zero.

The first of these radiative-acoustic modes is, however, somewhat mysterious, as it moves upstream through the flow at a speed faster than the sound speed in the gas. Owocki et al. (1986) demonstrated that Abbott’s approach is justified only for purely wave-like modes, in which the structure of the entire perturbation can be analytically continued from any finite part of it. This conclusion remains somewhat controversial, although it is analogous to recently observed laboratory phenomena such as the (apparently) superluminal propagation of signals in samples of Bose-Einstein condensates \[Wang et al. 2000\]. While information cannot propagate upstream because of the purely downstream direction of photon propagation, that does not prevent the phase and group velocities of wave-like perturbations being directed upstream in an unstable flow.

Using the downstream velocity gradient in the calculation of the Sobolev line force, however, allows information to travel upstream through the grid, and as a result genuine modes with the properties of Abbott waves arise in our numerical simulation. The speed of information propagation in these artificial waves is limited to $`v<\mathrm{\Delta }x/\mathrm{\Delta }t`$, which may be large if the time step is small on a coarse grid. The upstream velocity of these modes no longer requires the flow to be unstable, and in the outcome it is not. This can be illustrated, for example, by observing that the Abbott-like modes can allow the solution at the photosphere to adjust to satisfy the conditions at the CAK critical point, and hence the flow overall to relax to the steady-state CAK solution.

This stabilising feature of Abbott waves should, however, not be dismissed out of hand as a numerical artifact, as scattered radiation can lead to some upstream propagation of information. The downstream differencing scheme can be thought of as combining the (physical) upstream differencing scheme with a diffusion term

$$g(\rho _{n+1}\rho _n)/\mathrm{\Delta }x\alpha (\rho _{n+1}2\rho _n+\rho _{n1})/(\mathrm{\Delta }x)^2,$$
(34)
where $`\alpha =\mathrm{\Delta }x`$ gives downstream differencing, and $`\alpha =\mathrm{\Delta }x/2`$ gives a centred gradient. This diffusion term can then be interpreted as a first-order approximation to the effect of photons which scatter inwards at large radii and subsequently interact again at smaller radii (the scattered photons which do not subsequently interact are adequately treated by the absorption terms).

Specifically, the component $`\rho _n\rho _{n1}`$ corresponds to the outward scattering force from photons which interact at radii beyond the cell of interest: the downstream difference is allowable for these photons as they subsequently move inwards, causing an inward driving proportional to $`\rho _{n+1}\rho _n`$. We note that in this simple model, flow stability only results when the inward force resulting from scattering entirely cancels the direct force. The line-drag effect may be thought of as the result of just such a stabilizing influence, although at a different level to that implicit in fully downsteam gradient evaluation.

### 4.2 OCR force law

Initially the mean flow of the wind in one dimension was simulated from the above model. Although the model as described, is physically incomplete it is useful to explore the most simple case of a line-driven wind, i.e. unperturbed flow. Unperturbed here means that there is no extra perturbation inserted by hand, however the flow does contain features that arise from its unstable nature in the optically thin limit. This means that the flow is especially susceptible in the early wind to instabilities.

As was shown previously for the CAK force, the choice of the direction of velocity gradient calculation can effect the stability of the flow. Since VH-1 stores variables as cell averages, the determination of cell interface values for the velocity for use in the calculation of the radiative acceleration requires the choice of an averaging method. Two methods were used, a mixture of interpolation and extrapolation from upstream information, and interpolation using downstream and upstream information,

$$\text{mixture}\{\begin{array}{cc}v_\mathrm{u}\hfill & =\frac{1}{2}(v_n+v_{n1})\hfill \\ v_\mathrm{d}\hfill & =\frac{3}{2}(v_n\frac{1}{2}v_{n1})\hfill \end{array}$$

$$\text{interpolation}\{\begin{array}{cc}v_\mathrm{u}\hfill & =\frac{1}{2}(v_n+v_{n1})\hfill \\ v_\mathrm{d}\hfill & =\frac{1}{2}(v_n+v_{n+1})\hfill \end{array}$$

respectively. The subscripts refer to upstream (u) and downstream (d) values. The mixture method is the most physically consistent since it will not allow information to propagate upstream as a result of the numerical method. The interpolation method, since it uses both upstream and downstream grid cells, does permit information propagation upstream.

#### 4.2.1 Unperturbed flow

The results of these two methods of calculating the interface velocity for a simulation of the structure of an unperturbed flow are presented in Figure 3. It is apparent that while both methods are unstable, this is true to a far greater extent for the mixture method.

The instabilities are similar in form to the self-excited waves found by OCR. As in this previous paper, there is no clear mechanism associated with their seeding, although it seems probable that fluctuations caused by variations in the time step are amplified by the line-driven flow instability. In subsequent papers (Owocki 1991; Owocki & Puls 1996), Owocki & his collaborators have adopted an approach to the calculation of the scattered radiation, the Smooth Source Function approach, which allows the wind to recover the CAK solution when it is not perturbed at the photosphere.

However, as argued by Owocki & Rybicki (1984), where the scattering is small the wind should be intrinsically unstable throughout its volume. This is situation is similar to that often found in the study of conventional turbulent flows, and various approaches have been taken to treat it in that context. In one approach, simulations are regularized at the smallest resolvable lengthscales by a numerical hyperviscosity, but an appropriate small-scale stochastic input should be included throughout the grid to model the effect of sub-grid instabilities. We here take an alternative approach, in the spirit of large eddy simulations, allowing truncation errors to lead to perturbations in the flow. The exponential growth of the instabilities and their saturation at non-linear amplitude should mean that the overall properties of the flow are not strongly dependent on the non-physical nature of their seeding.

While the dependence of updated flow values on the grid data is more physically consistent in the mixture method, this method is also highly unstable. Some smoothing of the results is in fact desirable, in order to ensure that structures continue to be resolved by the numerical grid (this may be compared to the small smoothing of shocks necessary to prevent post-shock ringing in Godunov-type schemes). In real stellar winds, the line-drag effect reduces the net effect of the driving felt by these small wavelength disturbances at the base of the wind.

The interpolation method will therefore be consistently used in the remainder of this paper, and the flow shown by the solid line in Figure 3 as the initial condition for all the remaining simulations. This has the additional benefit that the development of finite perturbations may be studied in the relatively smooth region of the solution between $`z=1`$ and $`2`$ without being confused by the effects of rapidly-growing small-scale instabilities.

##### Statistical characterisation of the flow

The overall structure of the flow variations may be characterized by the statistical properties defined in equation (25). The results for the above model, are shown in Fig. 4a-c.

As one might expect, the velocity dispersion, $`v_{\mathrm{disp}}`$, is close to zero at the base of the wind and the clumping factor, $`f_{\mathrm{cl}}`$, is approximately unity. In this region there is little contribution to the flow structure from the noise excited waves, so the wind is smooth and has an almost constant velocity.

A minor peak in $`f_{\mathrm{cl}}`$ is reached at $`z=2`$, the radius where the noise waves begin to dominate the flow. Upstream of this point there are no dense shells but in the downstream flow the dense shells steepen into shocks. The flow contains many dense shells but it is dominated by the accompanying rarefactions. The noise in $`f_{\mathrm{cl}}`$ beyind $`z3`$ results from the extreme variability of the flow in this regime, so an average over $`100`$ time units has still not adequately sampled the behaviour of the flow in this region. Nevertheless, it is clear that the amount of clumping and velocity variation is increasing rapidly in the region.

Examination of $`v_{\mathrm{disp}}`$ in Fig. 4b shows that it too reaches a minor peak at $`z=1`$. Because this marks the point where the diffuse clumps form into dense shells bounded by forward and reverse shocks, $`v_{\mathrm{disp}}`$ decreases as the smooth rarefied region begins to dominate. This downturn is short lived as the rarefied regions broaden and the shocked gas moves downstream.

The behaviour of $`C`$ is a little surprising but not totally unexpected. A flow dominated by reverse shocks will have $`C1`$, since the velocity and density are anti-correlated. Fig. 4c shows this to be the case for the majority of the wind studied here. The initial portion of the wind moves sharply from correlated at the inner boundary to anti-correlated after a short distance, reaching a minimum at $`z=1.2`$. This implies that reverse shocks are the dominant structures, but it is unwise to ignore the forward shocks. Indeed, in the later stages of the flow each dense shell is bounded by both an upstream reverse shock and a downstream forward shock. These forward shocks, although small in amplitude, could begin to dominate the flow further downstream.

Comparing these results with Runacres & Owocki (2002), shows good agreement in all of the parameters except $`C`$. They observe a downstream value of $`C0`$ indicating that both forward and reverse shocks are equally prevalent. In the model here the reverse shocks control the flow to a greater extent than the forward shocks, borne out by the the slow rise, away from $`C=1`$ to $`C=0.5`$ at the edge of the grid.

##### Absorption spectrum

An indication of the line absorption spectrum which would be found from the flow can be found from the column density function, $`\eta (x)`$, at the edge of the grid. The high velocity part of the profile is highly variable, so we will first look at the time average of the structure.

In Fig. 5, we show this for the above model. Absorption profiles can be derived using

$$F(x)=\mathrm{exp}[\kappa _\mathrm{l}\eta (x)]$$
(35)
where $`\kappa _\mathrm{l}`$ is the opacity in a particular resonance line. However, as our models are only of a restricted region of the wind, such profiles would not be comparable in detail with observed P Cygni profiles.

The results can nevertheless be compared to the theoretical predictions of Abbott (1982). Looking at the plots for $`T_{\mathrm{eff}}=40,000\mathrm{K}`$ and $`50,000\mathrm{K}`$ in this paper, they show similarities to Fig. 5, most notably between $`v=0\text{}7`$. The sudden drop in flux at $`v=1.5`$ where the steep decrease in density stops and its behaviour becomes close to linear, representing the only high density, low velocity region in the flow. The smooth rise in Fig. 5 is due to the smoothly increasing velocity and decreasing density in the early part of the wind before noise begins to dominate. When a perturbation is applied to the base of the wind the smooth region of $`\eta `$ in our model is disturbed and becomes noisy showing the effect of perturbations on the spectrum.

The variation of $`\eta (x)`$ may be displayed in a trail plot, Fig. 6. The development of dense shells in velocity space is particularly clear in this plot. At low velocities, the trails are steep, with nearly constant gradient (i.e., the shells have a low, constant acceleration). As the shells reach a region where the radiation field is no longer strongly absorbed by underlying gas, the acceleration rapidly increses to a higher value, seen from the shallower gradients of the trails. Broad but gradually narrowing regions of relatively smooth opacity appear at lower velocities than some of the unresolved structures, corresponding to gas in the outer part of these shells ‘surfing’ on the enhanced radiation to the red side of the saturated absorption band. More complex structures develop, as shells interact.

Care must be taken when interpreting information from such a trail plot, since flow structures with the same velocity can appear at different places in the wind. This can make the origin of each trail ambiguous. Nevertheless, a trail that remains well defined for a significant time is the result of a single optically thick feature, since since the combination of differently spaced features would not have the same stability.

The points marked I, II, III have been identified with dense shells on the density plot Fig. 7, which is taken at $`t=0.414`$. In this case the three marked trails are the dense shells which stand out against the flow as being the most optically thick. It is clear from Figs. 6 and 7 that although the shell marked I is not the densest shell, it has the highest optical depth of the shells at this time. This is evident from its steep evolutionary path, indicating that the region accelerates slowly and as it recieves less driving.

The point marked as III in Fig. 6 indicates where the evolution of a dense lump changes from optically thin evolution to optically thick. The difference in the evolution is visible in two ways; the dark, steeper trail shows a marked difference to the earlier light, shallow trail of the optically thin evolution. This transition occurs because of minor shells colliding with the perturbation shell. As faster shells merge with the perturbation from upstream, the perturbation sweeps up slower shells downstream of it. These interactions increase the density of the shell and alter the optical properties making it optically thicker. As this behaviour continues the shell slows and at $`t=0.485`$ it merges with another shell and its once distinctive optical depth becomes muddled with that of a shell further upstream which is forming out of the noise. The perturbation shell does not disappear but from a spectral perspective it is indistinguishable from other parts of the flow, whilst remaining the densest shell in the region of flow under consideration.

#### 4.2.2 Optically thin perturbation

A perturbation is applied to the base of the wind which is 10 grid elements thick ($`\mathrm{\Delta }z=0.0375`$). The amplitude of the perturbation is a $`v_\mathrm{p}=10\times v_0`$ in velocity, and the change in density is $`\rho _\mathrm{p}=\rho _0/10`$, the subscript $`0`$ refers to background flow values; within this region, the velocity was increased to 10 times its previous value, and the density reduced by the same amount factor. These values were chosen to represent a typical short-wavelength ($`1/k\mathrm{}_{\mathrm{Sob}}`$) perturbation, to compare with the analytical work from Owocki & Rybicki (1984), discussed above.

In this limit the evolution of the perturbation should follow equation (24). The gas in the perturbation experiences a constant acceleration, $`\tau _0\delta v`$. In Fig. 8 the behaviour of the perturbation appears to be stable, as the gas in the perturbation is swept up into a dense shell which advects slowly across the grid. Behind this shell, the constant velocity gradient is characteristic of a free rarefaction.

Immediately ahead of the shell, the flow remains smoother than was the case in the unperturbed flow. It seems that the strong shock at the outer edge of the dense reduces the fetch available to the rapidly moving unstable outward radiative-acoustic waves, and so limits their amplitude in this region.

As the velocity gradient in the rarefaction decreases, instabilities form at the edges of the rarefaction, as can be seen in Fig. 9. Behind the dense shell, subsidiary shells form, while the smooth expansion of the gas upstream of the edge of the rarefaction also becomes highly disturbed. The optically thin subsidiary shells gather momentum from line driving and run up the diffuse wake of perturbation until they collide with the shock at the leading edge of the dense shell. This merging of shells is similar to that which occurs in the outer wind.

The dense shell is not subject to the formation of instabilities and is gently advected across the grid. This supports the asymptotic form, equation (23) which implies that this feature is stable for small wavelength perturbations: any change in the radiation force will only affect the phase of the perturbation, or the derivative of the velocity. The coupling of the radiation field with the dense shell is very weak which is evident in its unchanging amplitude, indicating that it is largely unaided by the driving force.

Shock structure begins to develop in the flow as an optically thin flow encounters a dense shell. A perfect example of a reverse shock can be seen in Fig. 10, a dense shell with a smoothly rarefied region preceding it. This figure also shows other shocks forming and interacting with the surrounding dense regions.

## 5 Numerical evolution of perturbations in two dimensions

The non-Sobolev approach described in the present paper can be naturally applied to multidimensional flows. In the present paper, we will concentrate on the simplest case, of the propagation of a single ray through plane-parallel flow. This is equivalent to studying a very narrow slice of the wind at a reasonable distance from the photosphere.

This has the obvious advantage that only a single ray calculation need be performed, parallel to one axis of the computational grid. It may be argued that this approach is little more than a series of one-dimensional simulations in separate but contiguous channels. However, while there is no tangential component to the radiation driving, we will see that the flows which are driven in this direction have a substantial effect on the structure of the wind.

The flows are studied in a similar way to the one-dimensional cases presented above. Perturbations are made to velocity field, added close to the base of the wind. A variety of widths and amplitudes of perturbation were modelled, but in the end the evolution of the flow always returns to structure of a similar form, so we will present here a small selection of typical results.

An initial two dimensional data set was obtained by filling each row of the simulation with identical data, derived from the unperturbed one dimensional run shown in Figure 3. Running the model without further perturbation, its development is similar initially to one-dimensional runs. However, as the simulation progresses these structures fragment in the tangential direction. It is to be presumed that these instabilities are seeded by small rounding errors in the underlying hydrodynamical scheme. Certainly the fragmentation initially appears nearly periodic.

The growth of these perturbations of the dense noise-excited shells may be attributed to Richtmeyer-Meshkov instability (a particular case of Rayleigh-Taylor instability which occurs when a shock interacts with a density discontinuity). As structure evolves from these instabilities, it has a very definite effect on the downstream flow. The result is a clumpy, highly structured flow as can be seen in Fig. 11.

If the gas is diffuse the instabilities do not evolve to macroscopic structure, because the strong coupling with the driving washes out the structure. Dense shells that have become optically thick, expand from gas pressure alone. This provides mechanisms by which both dense and diffuse shells can expand. Shocks are formed where the edges two such expanding shells collide, for example at (ii) in Fig. 12b. Depending on the expansion velocities of the shells involved, hydrodynamical instabilities may develop at the interface. The development of instabilities in the flow can fragment the resulting shell into multiply bullets and diffuse bubbles.

To study the development of these structures in more detail, we will now investigate the development of flows in which an explicit perturbation is included in the initial conditions.

### 5.1 Velocity perturbed flow

As a first example of the development of a perturbed two dimensional flow, we set the flow velocity to be $`20`$ per cent greater than its value at the equivalent depth in the surrounding flow in a region of $`10\times 40`$ grid elements (at a resolution of $`1`$ element $`=0.00375`$) at the base of the wind. Initially the disturbance expands smoothly through the background, and no significant effects appear on large scales. The matter in the initially perturbed region collects into a dense clump, with rarefied regions immediately upstream and downstream of it. The downstream rarefaction is simply the result of the increased velocity of the gas in the perturbed region, while the upstream perturbation is caused by the shadow of the perturbation. In this shadow, the gas experiences a decreased radiation force and so slows and is swept up by the clump.

As the perturbation increases in size, the changes of shadowing from gas shadowed by the perturbation (etc.) begin to have effects on the dense shells far downstream. The patch of the dense shell marked (i) in Fig. 12b, falls in the shadow of the perturbation, decelerates and becomes more diffuse, causing the first break in a shell structure. At this time an instability develops in the unshadowed parts of the dense shell where the optically thin faster flowing gas meets an optically thick slower dense shock. The disruption of the shells is dependent on whether they continue to contract as they evolve, which is in turn dependent on whether they are compressed by a shock behind or retarded by a dense region in ahead. If the shell expands slower than the perturbation evolves, the RM instability within the shell structure develops. If not the unstable structure diffuses away as the gas expands. This is the case in shell (ii) of Fig. 12c where the increasing radial width and decreasing density of the shell damps the instability, so that before it becomes significant the shell disperses. The shadowed region is unaffected by this expansion due to the alterations in spectrum that govern its motion. The apex of the shadowed region in the same plot forms a transient dense feature labelled (iii) which is created by the fast highly rarefied region immediately upstream of it, from which it accretes gas. As the flow evolves further, this clump disperses.

It is interesting to note that when part of a dense shell falls into the rarefied slipstream of the perturbation it experiences very little driving and is swept up by the dense clump of the perturbation. There is steady accretion from such incidents, increasing the density and extent of the perturbation. Were it not for this accretion the perturbation would become highly rarefied and disappear, but as is shown in Fig. 12d, part of the dense shell is accelerated and compressed into a high density bullet which rams through the perturbation.

Further downstream, Fig. 12e this feature develops into a bow shock/arrow shape and added to further by faster flowing gas accreting onto it from behind. This structure is preserved until the edge of the grid, but it is unlikely to remain in this form in the far downstream flow. When the wind has reached the terminal velocity adiabatic expansion will cause it to diffuse away.

Persistent features appear in the inner wind at the position of the initial perturbation. Instabilities seeded by these disrupt the formation of dense shells in the downstream flow. This is eventually smoothed out as the features spread in the tangential direction, Fig. 12f. However, as the flow is everywhere supersonic, the perturbations ought to have been advected downstream long before.

In an Eulerian hydrodynamical scheme, perturbations to smooth flow in a particular cell can only be removed exponentially with time (in essence by dilution). A sufficiently strong instability can decrease the exponential damping rate, or even force the perturbations to grow. This may be viewed as a limitation to the accuracy of these simulations – as indeed it would be a limitation to the accuracy of any Eulerian simulation with finite resolution (Lagrangian simulations would suffer from complementary problems as a result of sampling noise). However, in reality the flow will continually be affected by small scale perturbations from the photospheric flow. At late times, as the dense but noisy shells evolve they become highly structured as the hydrodynamics and the line-driving struggles for dominance in these regions. This flow is similar to that for unperturbed initial conditions, which suggests that it is not sensitive to the precise nature of the upstream perturbations, so long as they are small – the exponential growth of the instability erases the small details of its seeding.

### 5.2 Density perturbed flow

Contrasting behaviour is observed if a density perturbation is applied to the background flow. The motivation for this choice of perturbation was interest in a macroscopic persistent feature in the outer portion of the wind. The phenomenon of Co-rotating Interaction Regions (CIRs) has been discussed previously (e.g. Dessart & Chesneau 2002) but as yet their formation is not well understood.

Initially, the perturbation expands in the form expected for such an overpressured region from pure hydrodynamics. A diffuse bubble appears, with dense shocked plates at its upstream and downstream edges. However, as the flow develops, these plates move apart and become more condensed, while the material between them is accelerated by line driving into the downstream clump, or is swept up by the approach of the upstream edge or clump. The direction which material will move and which clump it will be swept into is dependent on whether the gas falls in the shadow of the upstream clump or not.

The downstream flow still shows evidence of the spectral interference from the perturbation. The background dense shells are disturbed initially on the line of the edges of the perturbation, since at these positions a range of velocities are present while in the body of the perturbation has a very narrow velocity spread; the downstream material is only disturbed if gas upstream of it blocks radiation in a relevant frequency range.

The subsequent evolution of the upstream and downstream components of the original perturbation is distinct from that found in the previous section, showing sensitive dependence of the details of the flow on the initial conditions. One obvious feature, shown in Fig. 13b, is that both clumps of the perturbation appear to be blown away by the outflowing wind whereas in Fig. 12b the perturbation forms bow shocks as it pushes through slower material ahead of it.

These results prompted an investigation of a higher density initial perturbation, so a density perturbation by a factor of 4 times was applied to the same initial conditions as above. As one might expect the evolution of this perturbation, as seen in Fig. 14a-f is similar to that seen in Fig. 13a-f. In this case the upstream edge of the expanding bubble withstands disruption by the flow and eventually forms a ’v’-shaped front, which after further evolution develops highly detailed structure in Fig. 14d as a shock passes through it leaving RM instabilities to develop.

The structure from the downstream edge of the bubble forms features in Fig. 14c-d that are reminiscent of the structure in Fig. 13c-d. The perturbation has a great effect on the persistent features which remain in the upstream flow, leading to the beginning of highly clumped flow in Fig. 14f.

### 5.3 Physical consequences

Discrete Absorption Components \[Kaper et al. 1999, DACs; e.g.\] are significant features seen in the absorption profiles of unsaturated lines. Current observation suggests that the kinematics of these features is dominated by a rotational component, but theory has not addressed the manner in which they form and their dynamics as they move through the wind. The results of the above sections suggest that features that are formed in the early parts of the wind, are dispersed rapidly if their propagation is slow enough to be destroyed by RM instabilities. DACs are formed in the early part of a real wind when a slowly rotating region suffers a compression from a neighbouring fast moving region. Higher density compressed portions of the wind are observed in line spectra as DACs.

The dense, long lived bullets which are produced in our simulations are a strikingly similar to those inferred in some empirical models of DACs \[Hamann et al. 2001, e.g.\]. Similar structures would be expected to form as a result of variations in the stellar flux due to dark surface features, which may naturally explain the origin of CIRs.

The ultimate aim of these models is to understand the nature of observed hot star winds. While the present results are by no means definitive, it is clear that cohesive absorbing knots are a natural feature of line-driven winds. However, further, more physically detailed, numerical simulations will be required before a detailed comparison to the observations of wind structures such as DACs can be performed.

## 6 Summary

We present models of the structure of unstable line-driven flows, calculated in two dimensions with a full hydrodynamical code and non-local radiation transport.

Our results compare well to those of more physically complete (but geometrically restricted) models in their common limits. The reproduction of the noise waves and the dominance of reverse shocks in the downstream structure is in good agreement with those seen by other non-Sobolev approaches (e.g. OCR, Runacres & Owocki 2002, Feldmeier 1995). This behaviour is borne out by the close similarity of the statistical descriptors of the flows we study to the results presented by Runacres & Owocki (2002). This suggests that the current model does indeed present a good description of the non-Sobolev behaviour of this wind.

The model presented here has limitations: 1) the gas is taken to be isothermal, which is reasonable at least in high density regions; 2) the spatial domain is restricted, effectively to a small region reasonably far from the photosphere; 3) the line profile function is taken to have a top-hat form rather than a Gaussian or Voigt profile; 4) the fiducial value of $`v_{\mathrm{th}}/c_\mathrm{s}=0.5`$ was used \[Owocki 1991, following\], although a rather lower value may in fact be more appropriate; 5) we have treated the radiation transfer in the pure-absorption approximation. We have argued, however, that none of these approximations will have a significant effect on the local structure of the line-driven winds, which is our principal interest in this paper: our work complements other studies in which these processes are treated in more physical detail in a more geometrically restricted domain.

In future work, we will include the line-drag effect of scattered radiation \[Lucy 1984\] \[Owocki & Puls 1996, formulated in the non-Sobolev method by\]. The asymmetry we have shown between forward and backward calculation of the Sobolev line force from the wind velocity gradient, could be of use in approximating the fore/aft asymmetry which produces line drag. Since this is a local effect, it may not be necessary to use a full non-Sobolev method for its calculation. Indeed a Sobolev approach may actually be more suitable than a global method for simulations with high degrees of imposed symmetry.

## Acknowledgments

This work was supported by PPARC through a research studentship (ELG) and an Advanced Fellowship (RJRW), and by NASA through grant NAG5-12020. We would like to thank Mark Runacres and Ian Stevens for several helpful discussions, and the referee, Stan Owocki, for a constructive and helpful report.