text:
# 1 Introduction

## 1 Introduction

The statistical theory of spectra  provides an appropriate method for examining the symmetry properties of quantum systems. Level statistics such as the nearest-neighbour spacing (NNS) distribution and the spectral rigidity $`\mathrm{\Delta }_3`$ are defined for a pure sequence of levels that have same quantum numbers. Usually, it is taken as granted that in a classically integrable system, which has as many integrals of motion (or quantum numbers) as the number of degrees of freedom, the levels are uncorrelated and so have a Poissonian NNS distribution. The pure level sequence in a time-reversal-invariant quantum system whose classical counterpart is chaotic, is successfully represented by a Gaussian orthogonal ensemble (GOE) of random matrices. A deviation of level statistics from the predictions of the GOE may be taken as an evidence that the sequence is not pure because of the existence of a hidden symmetry. The statistical theory of spectra has also been extensively used in the investigation of symmetry breaking because the destruction of a quantum number has a dramatic effect on the spectral fluctuations. Impressive are the studies of isospin mixing in light nuclei . Shriner et al. measured a ”complete” spectrum of low-lying states of 26Al involving states with isospins $`T=0`$ and 1 throughout the energy range covered by the data. These data offered a testing ground for studying the influence of isospin-symmetry breaking on the fluctuation properties of energy spectra (see, e.g., ). However, the limited data available in 26Al (142 levels) precluded a definite conclusion. In a recent experiment with monocrystalline quartz blocks, Ellegaard et al.  measured about 1400 well-resolved acoustic resonances. The properties of quartz allowed these authors to measure the gradual break-up of a point-group symmetry, which is statistically fully equivalent to the breaking of a quantum number like isospin. The spectra are obtained by externally tuning the symmetry breaking  , allowing a statistically significant investigation of the whole transition.

In the present paper we consider the effect of gradual symmetry breaking on the fluctuation properties of energy levels. When the symmetry is full the spectrum is divided into independent sequences, each corresponding to one of the eigenvalues of the symmetry operator. We consider two possible scenarios for breaking the symmetry. The first is provided by the Guhr-Weidenmüller model , which assumes that the perturbation acts equally on all states and that the matrix elements between states of different symmetry eigenvalues are random numbers of equal variances. In the second scenario, the perturbation acts differently on the unperturbed states and the symmetry breaking interaction mixes a limited number of the degenerate eigenstates that belong to different symmetry representations, and this number grows as we increase symmetry breaking. We are thus creating a new sequence of levels corresponding to states with no symmetry at the expense of the initial sequences. Consequently, as the fractional density of the no-symmetry sequence increases, the $`\mathrm{\Delta }_3`$ statistics of the total spectrum first increases, moving in the direction of the prediction of the Poisson statistics until the densities of all the sequences become equal, and then decreases to reach asymptotically the value given by the GOE. Section 3 shows that the behaviour described by the second scenario is indeed present in the data of Ellegaard et al. . The summary and conclusions of this work are given in Sect. 4.

## 2 Symmetry breaking in a chaotic system

We consider a chaotic system that conserves a given symmetry. In such a system, the Hamiltonian assumes a block-diagonal form :

$$H_S=\mathrm{diag}(H_1,H_2,\mathrm{},H_n),$$
(1)
where each sub-block corresponds to states of certain quantum numbers. Thus, when the symmetry is conserved, the levels are divided into a number of ”pure” sequences, each described by a pure GOE. However, the statistical properties of the total spectrum that combines these sequences are no longer of the GOE type. In the limit of infinite number of sequences combined, the spectrum obeys the Poisson statistics. For a finite number of sequences, the fluctuation properties are intermediate between the Poisson and GOE statistics. In the following, we shall restrict our consideration to the case when the symmetry under investigation is represented by an operator that has two possible eigenvalues, such as parity, or isospin 0 and 1 as considered in Refs .

The NNS distribution of a spectrum resulting from a random superposition of $`n`$ independent sequences is calculated, e.g., by Berry and Robnik  and in Mehta’s book . If the level density of the $`i`$th sequence is $`\rho _i`$, and if the NNS distribution of levels of this sequence is $`P_i(x_i)`$, where $`x_i=f_is`$, $`f_i=\rho _i/\mathrm{\Sigma }\rho _i,`$ and $`s`$ is the NNS normalised to a unit mean, then the NNS distribution of the mixed sequence is given by

$$P(s)=E(s)\left\{\underset{i}{}f_i^2\frac{P_i\left(f_is\right)}{E_i\left(f_is\right)}+\left[\underset{i}{}f_i\frac{1W_i(f_is)}{E_i(f_is)}\right]^2\underset{i}{}\left[f_i\frac{1W_i(f_is)}{E_i(f_is)}\right]^2\right\},$$
(2)
where $`E(s)=_{i=1}^nE_i(f_is),`$ $`E_i(x_i)=_{x_i}^{\mathrm{}}\left[1W_i\left(x\right)\right]𝑑x`$, and $`W_i\left(x_i\right)=_0^{x_i}P_i(x)𝑑x.`$ In particular, if all the $`n`$ individual sequences have the same level densities, so that $`f_i=1/n`$, and if the NNS distribution in each is a Wigner distribution

$$P_i(x_i)=\frac{\pi }{2}x_ie^{\frac{\pi }{4}x_i^2}$$
(3)
which is a good approximation for the NNS distribution of a GOE, then (2) becomes

$$P_n(s)=\frac{1}{n}\left[\mathrm{erfc}\left(\frac{s\sqrt{\pi }}{2n}\right)\right]^nQ_n(s)\left[\frac{\pi }{2n}s+(n1)Q_n(s)\right],$$
(4)
where $`Q_n(s)=\mathrm{exp}(\pi ^2s^2/4n^2)/`$erfc$`\left(s\sqrt{\pi }/2n\right)`$ and erfc$`(x)`$ is the complementary error function. This distribution has a shape intermediate between those of the Wigner and the Poisson distributions. The Brody distribution ,

$$P_\mathrm{B}\left(s\right)=c_\beta s^\beta \mathrm{exp}\left(\frac{c_\beta }{\beta +1}s^{\beta +1}\right)\mathrm{with}\mathrm{c}_\beta =\frac{\mathrm{\Gamma }^{\beta +1}\left(1/\left(\beta +1\right)\right)}{\beta +1},$$
(5)
accurately reproduces the spectra of many two-dimensional systems with mixed regular-chaotic dynamics. As we vary the parameter $`\beta `$ from 0 to 1 the statistics changes from Poissonian to GOE.

The NNS distribution contains information about the spectrum in a short range, not exceeding a maximum of three mean level spacings. Long range information is provided by various higher-order correlation functions . The most popular among these is the spectral rigidity $`\mathrm{\Delta }_3`$ of Dyson and Mehta . Semiclassical arguments  show that, for a generic integrable system whose spectrum satisfies Poisson statistics

$$\mathrm{\Delta }_{3,\mathrm{Poisson}}(L)=\frac{L}{15}\mathrm{for}\mathrm{L}\mathrm{L}_{\mathrm{max}},$$
(6)
where $`L_{\mathrm{max}}=2\pi \mathrm{}/DT_{\mathrm{min}}`$, $`D`$ is the mean level spacing and $`T_{\mathrm{min}}`$ is the period of the shortest classical orbit. For time-reversal-invariant chaotic system the semiclassical theory gives the following asymptotic expression

$$\mathrm{\Delta }_{3,\mathrm{GOE}}(L)=\frac{1}{\pi ^2}\mathrm{ln}L0.007\mathrm{for}1\mathrm{L}\mathrm{L}_{\mathrm{max}}.$$
(7)
In all cases $`\mathrm{\Delta }_3`$ saturates to a non-universal value at $`LL_{\mathrm{max}}`$. The random matrix theory obtains an analytical expression for $`\mathrm{\Delta }_{3,\mathrm{GOE}}(L)`$ that involves double integration . We find it more suitable to parametrise this expression for a GOE, guided by Eq. (7), in the form

$$\mathrm{\Delta }_{3,\mathrm{GOE}}(L)=\frac{1}{\pi ^2}\left(1e^{aL}\right)\left[\mathrm{ln}(L+b)+c\right],$$
(8)
with $`a=8.241,b=0.944`$ and $`c=0.06724`$ being the best fit parameters for values of $`L`$ in the range of $`0L50`$. Seligman and Verbaarschot  suggested the following expression for the spectral rigidity of a spectrum resulting from a random superposition of $`n`$ independent sequences with fractional level densities $`f_i`$

$$\mathrm{\Delta }_3(L)=\underset{i=1}{\overset{n}{}}\left[\mathrm{\Delta }_{3,i}(f_iL)\right],$$
(9)
where $`\mathrm{\Delta }_{3,i}(L)`$ is the spectral rigidity for the $`i`$th sequence. This proposal is not surprising because $`\mathrm{\Delta }_3(L)`$ is essentially a variance.

We now consider symmetry breaking in the chaotic system using the Guhr-Weidenmüller model . The model suggests that the effect of the symmetry-breaking perturbation on all the unperturbed states is equal on the average. This suggests defining the Hamiltonian in the form:

$$H=\left[\begin{array}{cc}H_0& 0\\ 0& H_0\end{array}\right]+\alpha \left[\begin{array}{cc}0& W\\ W^T& 0\end{array}\right],$$
(10)
where $`H_0`$ represents a GOE and $`W`$ is a random matrix whose elements have Gaussian distributions with zero means and variances equal to those of the non-diagonal elements of $`H_0`$. When $`\alpha =0,`$ the symmetry is conserved and the Hamiltonian has the block diagonal form of Eq. (1) with $`n=2`$. The spectrum consists of an independent superposition of two level sequences having equal density. The NNS distribution $`P(s)`$ is given by Eq. (4) and, in particular, $`P(0)=1/2.`$ The spectral rigidity is given by (9) with $`f_1=f_2=1/2`$. When we allow $`\alpha `$ to increase, the system gradually evolves toward a single GOE. To show how this transition proceeds, we considered an ensemble of ten 200$`\times 200`$ Hamiltonian matrices of the form (10) for a fixed value of the perturbation strength $`\alpha `$. We numerically diagonalised these matrices. We obtained the NNS distribution and spectral rigidity for each of the resulting spectra, and measured their deviation from the GOE statistics by comparing with the Brody distribution, Eq. (5), for the NNS distribution and with the Seligman and Verbaarschot expression, Eq. (9), for the $`\mathrm{\Delta }_3`$ statistics where $`\mathrm{\Delta }_{3,1}`$ is given by the Poisson distribution (6) and $`\mathrm{\Delta }_{3,2}`$ by the GOE (8), so that

$$\mathrm{\Delta }_{3,\mathrm{SV}}(L)=q_{\mathrm{SV}}\frac{L}{15}+\mathrm{\Delta }_{3,\mathrm{GOE}}\left[(1q_{\mathrm{SV}})L\right],$$
(11)
We then repeated the procedure for other value of $`\alpha `$. The best fit of the tuning parameters $`\beta `$ and $`q_{SV}`$ obtained for $`P(s)`$ and $`\mathrm{\Delta }_3(L)`$ corresponding to different values of the perturbation strength $`\alpha `$ are shown in Fig. 1. The error bars measure the changes of the parameters that increase the minimum $`\chi ^2`$ values by a factor of 2. The monotonous nature of the curves clearly shows that the transition from symmetry to no-symmetry in the Guhr-Weidenmüller model is gradual and one-directional.

The opposite extreme of the Guhr-Weidenmüller model is the case when the eigenfunctions will become either very weakly or very strongly perturbed. For those states in which the perturbation is weak enough, the wave function can in principle be characterised by one of the two symmetry eigenvalues. These will be the states for which the splitting of two degenerate levels belonging to the two symmetry representations can be neglected. The corresponding energy levels will again form two independent sequences with fractional level densities $`f_1=f_2=f/2,`$ while the levels of the strongly perturbed eigenstates will form a new sequence with density $`f_3=1f`$ which will grow at the expense of the other two as the symmetry-breaking perturbation increases. The NNS distribution will then be given, using (2) and (3), by

$`P(s)`$ $`=`$ $`{\displaystyle \frac{\pi }{8}}f^3se^{\pi f^2s^2/16}\mathrm{erfc}\left({\displaystyle \frac{\sqrt{\pi }}{4}}fs\right)\mathrm{erfc}\left[{\displaystyle \frac{\sqrt{\pi }}{2}}\left(1f\right)s\right]`$ (12)
$`+{\displaystyle \frac{1}{2}}f^2e^{\pi f^2s^2/8}\mathrm{erfc}\left[{\displaystyle \frac{\sqrt{\pi }}{2}}\left(1f\right)s\right]`$
$`+{\displaystyle \frac{\pi }{2}}\left(1f\right)^3se^{\pi (1f)^2s^2/4}\left[\mathrm{erfc}\left({\displaystyle \frac{\sqrt{\pi }}{4}}fs\right)\right]^2`$
$`+2f(1f)e^{\pi \left(5f^28f+4\right)s^2/16}\mathrm{erfc}\left({\displaystyle \frac{\sqrt{\pi }}{4}}fs\right).`$
for spacings greater than a certain lower limit. If we ignore the latter condition, the extrapolation of $`P(s)`$ into the origin yields

$$P(0)=f\left(2\frac{3}{2}f\right),$$
(13)
We note that $`P(0)`$ increases with decreasing $`f`$ from a value of 1/2 at $`f=1`$, reaches a maximum value of 2/3 at $`f=2/3`$, and then decreases monotonously until it vanishes at $`f=0`$. Equations (12) and (13) show that the evolution of the shape of the NNS distributions with increasing symmetry breaking, and decreasing $`f`$, is not a straightforward transition from the 2-GOE behaviour to that of a GOE. Indeed, by decreasing the value of $`f`$ we are removing some of the levels from the initial two symmetry-invariant sequences to form the third sequence that represents the states without symmetry. The distribution will look more regular when we decrease $`f`$ until the level densities of the three sequences become equal. If we further increase the symmetry breaking, the initial sequences become thinner until the whole spectrum forms a single GOE sequence.

We can follow this transition better by calculating the spectral rigidity which is, for this system, given by

$$\mathrm{\Delta }_3\left(L\right)=2\mathrm{\Delta }_{3,\mathrm{GOE}}\left(\frac{fL}{2}\right)+\mathrm{\Delta }_{3,\mathrm{GOE}}\left[\left(1f\right)L\right].$$
(14)
In Fig. 2, the spectral rigidities corresponding to $`f=1`$, 0.9, 0.667, 0.333, 0.1 and 0 are given by the curves labelled by A, B, C, D, E and F, respectively, while the curve P corresponds to a regular spectrum. The figure clearly shows that $`\mathrm{\Delta }_3\left(L\right)`$ increases as symmetry breaking increases, and thus evolves towards the regular shape until $`f`$ reaches the value of 2/3. Only then starts $`\mathrm{\Delta }_3\left(L\right)`$ a route leading to the GOE curve.

## 3 Comparison with acoustic resonance spectra

In this section, we show that the proposed application of the Berry-Robnik theory to symmetry breaking leads to a satisfactory description of the acoustic resonance spectra of monocrystalline quartz blocks, measured by Ellegaard et al. . Crystalline quartz exhibits $`D_3`$ point-group symmetry about the crystal’s $`Z`$ (optical) axis, and three two-fold rotation symmetries about the three $`X`$ (piezoelectric) axes; the latter three axes lie in a plane orthogonal to the $`Z`$ axis and sustain angles of 120 degrees with respect to each other. Ellegaard et al. used rectangular blocks of dimensions $`14\times 25\times 40`$ mm$`^3,`$ cut in such a way that all symmetries are fully broken except a two-fold ”flip” symmetry about one of the three $`X`$ axes. A gradual breaking of this symmetry is then achieved by removing successively bigger octants of a sphere from one corner, thereby creating a three-dimensional Sinai billiard. The acoustic spectrum is measured for 8 radii of the removed octant, providing data of high statistical significance for the study of the gradual breaking of a point-group symmetry.

Figure 3 shows the spectral rigidity for the cases when the symmetry is present (open circles), and when it is violated by removing a tiny octant of radius $`r=0.5`$ mm (closed circles). The straight line labelled P is for an integrable system. We note that data for the case of tiny symmetry violation are systematically higher than the ones for the case of conserved symmetry. This agrees with the prediction of the second scenario proposed above, in which the system that undergoes a gradual symmetry breaking looks as if it has become more regular in the early stages of the transition.

Ellegaard et al.  concluded from the rise of $`\mathrm{\Delta }_3`$ over the theoretical prediction for a superposition of two GOE sequences that the quartz block with conserved flip symmetry has much in common with a pseudo-integrable system. They supported this conclusion by means of other measurement. Pseudo-integrable (PI) systems are non-integrable systems, yet non-chaotic. Two-dimensional PI systems, exemplified by particles moving in a planar polygonal enclosure with rational angles , share with the integrable system that the motion of trajectories in the phase space is restricted to two-dimensional compact surfaces. But these invariant surfaces are multi-connected and not tori as in the case of integrable systems. Numerical studies of the $`\pi /3`$ rhombus billiard  show that the spectral fluctuation properties are intermediate between those of a Poisson ensemble and a GOE. Bogomolny et al.  studied the level statistics of PI systems using a modified version of Dyson’s stochastic Coulomb gas model , in which the interaction between particles is restricted only to the nearest neighbours. The NNS distribution obtained in  is given by

$$P_{\mathrm{PI}}(s)=4se^{2s}.$$
(15)
This distribution is in agreement with the numerical calculation for the rational billiards . Unfortunately, the formula obtained by these authors for the two-point correlation function

$$R_2(s)=14e^{4s}$$
(16)
leads to an expression for the $`\mathrm{\Delta }_3`$ statistics that asymptotically increases as $`L/30`$, which is not consistent with the results of the numerical experiments. Following several authors, e.g. , we shall represent the $`\mathrm{\Delta }_3`$ statistics by Seligman-Verbaarschot  formula (9), obtained for a superposition of Poisson and GOE level sequences,

$$\mathrm{\Delta }_{3,\mathrm{PI}}(L)=q\frac{L}{15}+\mathrm{\Delta }_{3,\mathrm{GOE}}\left[(1q)L\right],$$
(17)
where $`q`$ is a fitting parameter. Biswas and Jain  found $`q`$ to be equal to 0.2 for the $`\pi /3`$ rhombus billiard.

We now apply the second scenario, which has been proposed in the previous section, to the symmetry breaking in a PI system having a symmetry represented by a quantum number that takes one of two possible values. When the symmetry is conserved, the spectrum is an independent superposition of two sequences corresponding to the two representations of the symmetry. The NNS distribution is accordingly represented by (2) and the $`\mathrm{\Delta }_3`$ statistics by (9) with $`n=2`$, $`P_i\left(x_i\right)`$ given by (15) and $`\mathrm{\Delta }_{3,i}(L)`$ by (17). We now again assume that as the symmetry violation starts, the eigenfunctions of the system will be divided into two classes. The first constitutes a fraction $`f`$ of the eigenfunctions in which the perturbation due to symmetry breaking can be ignored. The levels of these eigenstates are divided into two independent sequences exactly as before switching the symmetry breaking on. The second is the class of the strongly perturbed states that have completely lost the symmetry. This latter class will here be modelled by a GOE. Therefore, the total spectrum will be composed of 3 independent sequences, two with PI statistics and each having a fractional density of $`f/2`$, and one with GOE statistics and fractional density $`1f`$. During the symmetry-breaking transition, the NNS distribution is given by

$`P(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}f^2\left(2f^2s^2+4fs+1\right)e^{2fs}\mathrm{erfc}\left[{\displaystyle \frac{\sqrt{\pi }}{2}}(1\mathrm{f})\mathrm{s}\right]`$ (18)
$`+{\displaystyle \frac{1}{8}}(1f)(2+fs)\times \left[\pi fs^2(1f)^2+2s\left\{\left(\pi +4\right)f^22\pi f+\pi \right\}+8f\right]`$
$`\times e^{\frac{\pi }{4}(1f)^2s^22fs},`$
which is obtained by substituting Eqs. (3) and (15) into Eq. (2). The $`\mathrm{\Delta }_3`$ is given by

$$\mathrm{\Delta }_3\left(L\right)=2\mathrm{\Delta }_{3,\mathrm{PI}}\left(\frac{fL}{2}\right)+\mathrm{\Delta }_{3,\mathrm{GOE}}\left[\left(1f\right)L\right].$$
(19)
We now compare our model with the acoustic resonance spectra measured by Ellegaard et al. . We start by estimating the parameter $`q`$ in the expression (17) for the $`\mathrm{\Delta }_3`$ statistics of the PI system. We do this by fitting the same statistics for the spectrum of the complete quartz cube to Eq. (19) with $`f=1`$ and obtain $`q=0.082`$. We then compare Eq. (19) with the spectral rigidities of the other spectra considering $`f`$ as a free parameter. The best-fit values of $`f`$ for spectra corresponding to different radii of the removed octants are given in Table 1. We use these values to calculate the NNS distributions. The results of the calculation are compared with the experimental data in Fig. 4. We see that the proposed model presents a satisfactory description of the whole transition by varying a single parameter. The only disagreement is between the calculated and measured values of $`P(s)`$ at small values of $`s`$ as we have already expected. The symmetry breaking decreases the probability of finding degenerate levels sharply leading to the observed dip at small $`s`$ in the spacing distributions . This dip is followed by an overshoot to restore normalisation. The width of this dip, which is about 1/10 of the mean level spacing, provides an estimate for the minimum level-splitting below which two levels can be regarded as approximately degenerate.

## 4 Summary and conclusions

Conventional models for symmetry breaking in chaotic systems assume that the non-diagonal matrix elements of the symmetry breaking perturbation are statistically equivalent. They predict a monotonous transition from a mixed statistics to the GOE as one switches on and increases the symmetry-breaking perturbation. This is not the behaviour of the acoustic resonance spectra of quartz measured by Ellegaard et al . We propose another scenario for symmetry breaking in which the levels are approximately separated into two groups, one with and one without symmetry. This happens when the expansion of any of the wavefunctions in terms of the eigenfunctions of the symmetry operator is either dominated by a single term or composed of statistically equivalent contributions from all its component. Considering a system in which a symmetry that has two representations is conserved, we represent the spectrum as a superposition of two independent sequences, one for each symmetry representation. Now, as we switch on the symmetry breaking perturbation, we remove a number of levels from these two level sequences into a third sequence. The latter corresponds to the strongly perturbed states in which the symmetry has virtually disappeared. Because we have increased the number of spectral partitions, the spectrum will look as if it became more regular by the small violation of the symmetry. The apparent regularity will continue until the three sequences become equally populated. The further increase of symmetry breaking will reverse the evolution of the shapes of the spectral statistics toward the GOE shapes. Namely this is the behaviour found by inspection of the acoustic resonance spectra of quartz measured by Ellegaard et al . We see, for instance, from Fig. 4 that it is possible to describe the evolution of the NNS distribution and the spectral rigidity during the transition from a fully conserved symmetry to a fully violated one by varying a single parameter, measuring the fractional density of states which are practically unaffected by symmetry violation.

Table 1. Best fit values for the fraction $`f`$ of states not affected as a result of violating the point-group symmetry of a quartz rectangular block by removing spherical octants of various radii $`r`$. The values of $`f`$ are evaluated by fitting the spectral rigidities $`\mathrm{\Delta }_3\left(L\right)`$ of Fig. 4. to Eq. (19).

| $`r`$ | 0 | 0.5 mm | 0.8 mm | 1.1 mm | 1.4 mm | 1.7 mm | huge defocusing |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $`f`$ | 1 | 0.49 | 0.37 | 0.28 | 0.23 | 0.14 | 0 |

Figure Caption

Figure 1. The monotonous dependence of (a) the parameter $`\beta `$ in the Brody distribution (Eq.(5)) for the NNS distributions $`P_\mathrm{B}(s)`$ and (b) the parameter $`q_{SV}`$ in the Seligman and Verbaarschot expression (Eq.(11)) for the spectral rigidity $`\mathrm{\Delta }_3\left(L\right)`$ that interpolate between the Poisson and GOE statistics, obtained by a numerical experiment on a chaotic system undergoing a gradual breaking of a symmetry according to the Guhr-Weidenmüller model .

Figure 2. The spectral rigidity $`\mathrm{\Delta }_3\left(L\right)`$ for a chaotic system undergoing a gradual breaking of a symmetry according to the scenario proposed in the text. The results of the calculation when the fraction of states remaining nearly symmetry invariant $`f=1`$, 0.9, 0.667, 0.333, 0.1 and 0 are given by the curves labelled by A, B, C, D, E and F, respectively, while the curve P corresponds to a regular spectrum.

Figure 3. The spectral rigidity $`\mathrm{\Delta }_3\left(L\right)`$ measured in  for the quartz block when the symmetry is present (open circles), and when it is violated by removing a tiny octant of radius $`r=0.5`$ mm (closed circles). The straight line labelled P is for an integrable system.

Figure 4. The NNS distributions $`P(s)`$ and spectral rigidities $`\mathrm{\Delta }_3\left(L\right)`$ for the different radii of the octant removed from a quartz block: (a) $`r=0`$, the flip symmetry is fully conserved, (b) $`r=0.5`$ mm, (c) $`r=0.8`$ mm, (d) $`r=1.1`$ mm, (e) $`r=1.4`$ mm, (f) $`r=1.7`$ mm, (x) the block with the huge defocusing structure. The data are from Ref. , while the curves are the theoretical curves using Eqs. (18) and (19). Values of the parameter $`f`$ that measures the fractional densities of states that keep the symmetry are given in Table 1.


format: text

text:
# Classical calculation of high-order harmonic generation of atomic and molecular gases in intense laser fields

## I Introduction

The development of the high-power femtosecond laser has stimulated the investigation of the multi-photon processes of atoms and molecules interacting with intense laser fields \[1-9\]. Recently, there are many theoretical and experimental references about these multiphoton process. Within a proper intensity region, lots of odd harmonics of laser are generated by atomic and molecular gases \[1-5,14\]. The harmonic structure distributes as a plateau which is cut off at a special high-order harmonic. When the intensity increases, the ionization channel is opened and the high-order harmonics disappear. The above threshold ionization (ATI) occurs when the laser intensity is sufficient power (above $`10^{13}W/cm^2`$). Various structures (plateau, angular distribution, etc.) in ATI spectra have been detailed in Ref.\[6-8\]. The above threshold dissociation (ATD) and dissociation-ionization of molecular systems are also reported \[9-10\]. Studying the laser-matter interaction deeply, not only can obtain new knowledge of the interacting mechanism, but also can provide widely application in the generation of high-order coherent harmonics, $`X`$-rays laser and $`\gamma `$-rays laser.

The classical dynamics of most laser-driven systems is generally chaotic, due to the existence of nonlinearity. Chaos usually manifests itself as some control parameters (initial energy, laser intensity, laser frequency, etc.) are varied. The microscopic systems, in particular those involving atoms and molecules, are governed by a Hamiltonian. To study these systems are of great importance in the context of quantum-classical correspondence. At the investigated high power laser intensity $`(10^{13}\mathrm{~}10^{15}W/cm^2)`$, the electric field of the laser is equal in strength to the Coulomb field of the nuclei , and the laser field can not be looked as a perturbation to the field-free system, for the states themselves are no longer independent of the laser. Generally, the microscopic systems (atoms or molecules) are intrinsically quantum mechanical systems, thus they must be described by quantum mechanics. However, due to the presence of intense laser, the exact calculation even numerical simulation based upon quantum mechanics tend to be very difficult to perform. Fortunately, classical approach to these similar problems is useful for providing physical insight into dynamics processes \[10-13\]. Classical chaos associating with the microwave ionization of atomic hydrogen reveals that the detailed mechanism of atomic ionization in terms of transport in phase space \[11-12\]. Classical prediction for scaled frequencies has been verified by quantum calculations \[12-13\], and, in turn, has confirmed the dynamical significance of classical chaos. The harmonic generation (HG) of laser driven hydrogen atoms is simulated with classical Monte-Carlo method , the numerical results are qualitatively consistent with the quantum mechanic results and the experimental observation. For the molecular systems, the classical calculation is a good first step, since length and energy scales are often large enough for classical mechanics to be at least approximately valid. Using softened potential model, the classical dynamics of the one-dimensional hydrogen molecular ion $`H_2^+`$ interacting with an intense laser pulse are detailed \[15-16\].

There are two important and realistic examples of laser driven systems, where nonlinear dynamics has played a major role, which have been traditionally studied. One is the laser driven hydrogen atom that is the fundamental system in atomic physics, the other is the laser driven hydrogen molecular ion $`H_2^+`$, which is the fundamental diatomic molecular system in molecular physics. Bellow, we shall show our classical calculation of the high-order harmonic generation of these two fundamental systems. The outline of this paper is as follows. The regularized model and the corresponding initial microcanonical distribution are presented in the next section. In section III, we show the numerical results in details. A briefly summary and discussion is given out in the last section.

## II Model and initial microcanonical distribution

### A Regularized model of laser driven hydrogen atom

In this article, we consider a classical hydrogen atom, with an infinite mass nucleus fixed at the origin of coordinates, interacting with a high intense laser field which is linearly polarized along the $`zaxis`$, and with the electric field component $`\epsilon (t)`$. Thus the Hamiltonian in atomic units in Cartesian coordinates is

$`H`$ $`=`$ $`H_0+H_{i,}`$ (1)
$`H_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}p^2{\displaystyle \frac{1}{r}},H_i=z\epsilon (t).`$ (2)
where, $`r=\sqrt{x^2+y^2+z^2}`$ and $`p^2=p_x^2+p_y^2+p_z^2`$. Apparently, the above Hamiltonian has a Coulomb singularity corresponding to electron-nucleus collision. To remedy the singularity, we introduce the parabolic coordinates $`(u,v,\varphi )`$

$$x=uv\mathrm{cos}\varphi ,y=uv\mathrm{sin}\varphi ,z=(v^2u^2)/2.$$
(3)
and a new fictive time scale $`\tau `$

$$dt/d\tau =G(u,v)=v^2+u^2.$$
(4)
In order to implement regularization, following the notation of Szebenhely , regarding the motion time $`t`$ and the negative total energy $`E`$ as generalized coordinate and generalized momentum respectively, then the Hamiltonian function in extended space becomes into

$$H^{}=HE0=H^{}(x,y,z,t;p_x,p_y,p_z,E)$$
(5)
For the sake of obtaining the regularized Hamiltonian, we introduce the third category generating function $`F_3(p_x,p_y,p_z,E;u,v,\varphi ,t)`$, then obtain

$$x(u,v,\varphi )=\frac{F_3}{p_x},y(u,v,\varphi )=\frac{F_3}{p_y},z(u,v,\varphi )=\frac{F_3}{p_z},t(u,v,\varphi )=\frac{F_3}{E}.$$
(6)
So the third category generating function $`F_3(p_x,p_y,p_z,E;u,v,\varphi ,t)`$ is in the form of

$`F_3`$ $`=`$ $`xp_xyp_yzp_z+Et,`$ (7)
$`=`$ $`uv\mathrm{cos}\varphi p_xuv\mathrm{sin}\varphi p_y{\displaystyle \frac{1}{2}}(v^2u^2)p_z+Et.`$ (8)
The new momenta $`(p_u,p_v,p_\varphi ,E)`$ are related with old momenta $`(p_x,p_y,p_z,E)`$ by

$`p_u`$ $`=`$ $`{\displaystyle \frac{F_3}{u}}=v\mathrm{cos}\varphi p_x+v\mathrm{sin}\varphi p_yup_z,`$ (9)
$`p_v`$ $`=`$ $`{\displaystyle \frac{F_3}{v}}=u\mathrm{cos}\varphi p_x+u\mathrm{sin}\varphi p_yvp_z,`$ (10)
$`p_\varphi `$ $`=`$ $`{\displaystyle \frac{F_3}{\varphi }}=uv\mathrm{sin}\varphi p_x+uv\mathrm{cos}\varphi p_y,`$ (11)
$`E`$ $`=`$ $`{\displaystyle \frac{F_3}{t}}=E.`$ (12)
Then the regularized Hamiltonian can be expressed as

$`K`$ $`=`$ $`{\displaystyle \frac{dt}{d\tau }}H^{}=G(u,v)(HE)0,`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}[p_u^2+p_v^2+(u^2+v^2)p_\varphi ^2]2E(u^2+v^2){\displaystyle \frac{1}{2}}(v^4u^4)\epsilon (t).`$ (14)
If we define

$`K_u`$ $`=`$ $`{\displaystyle \frac{1}{2}}[p_u^2+u^2p_\varphi ^2]1Eu^2+{\displaystyle \frac{1}{2}}u^4\epsilon (t),`$ (15)
$`K_v`$ $`=`$ $`{\displaystyle \frac{1}{2}}[p_v^2+v^2p_\varphi ^2]1Ev^2{\displaystyle \frac{1}{2}}v^4\epsilon (t).`$ (16)
In the field-free case, $`\epsilon (t)=0`$, $`K_u(=K_v)`$ is the component of the Laplace-Runge-Lenz vector along $`zaxis`$, i.e., $`K_u=K_v=R_z`$. The equations of motion can be derived from the above regularized Hamiltonian, as the following:

$$\frac{d\varphi }{d\tau }=\frac{K}{p_\varphi }=p_\varphi (u^2+v^2),\text{ }\frac{dp_\varphi }{d\tau }=\frac{K}{\varphi }=0,$$
(17)
$$\frac{du}{d\tau }=\frac{K}{p_u}=p_u,\text{ }\frac{dp_u}{d\tau }=\frac{K}{u}=2Eu+p_\varphi ^2u^32u^3\epsilon (t),$$
(18)
$$\frac{dv}{d\tau }=\frac{K}{p_v}=p_v,\text{ }\frac{dp_v}{d\tau }=\frac{K}{v}=2Ev+p_\varphi ^2v^3+2v^3\epsilon (t),$$
(19)
and

$$\frac{dE}{d\tau }=\frac{K}{t}=\frac{1}{2}(u^4v^4)\frac{d\epsilon (t)}{dt},\text{ }\frac{dt}{d\tau }=\frac{K}{E}=G(u,v).$$
(20)
Apparently, $`p_\varphi `$ is a constant which corresponding to the component of angular momentum along $`zaxis`$.

### B Regularized model of laser driven hydrogen molecular ion $`H_2^+`$

Within Born-Oppenheimer approximation, we can assume that two protons are fixed at the positions $`A`$ and $`B`$, with a distance $`R`$ away from each other, in the locations $`(0,0,R/2)`$ and $`(0,0,R/2)`$ of the Cartesian coordinates system. A single electron at $`(x,y,z)`$ is subject to the Coulomb attraction of both protons and the interaction of the laser. Let $`r_{A\text{ }}`$ be its distance from $`A`$, and $`r_{B\text{ }}`$ be its distance from $`B`$. Within the dipole approximation, the Hamiltonian in atomic units is

$`H`$ $`=`$ $`H_0+H_{i,}`$ (21)
$`H_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}p^21/r_A1/r_B,H_i=z\epsilon (t).`$ (22)
where, $`\epsilon (t)`$ is the electric field of the laser pulse, $`H_i`$ is the interacting Hamiltonian. Apparently, the above Hamiltonian has singularities at points $`r_{A\text{ }}=0`$ and $`r_B=0`$, which corresponds to electron-proton collision. To overcome this barrier, regularization has to be performed. Define the new coordinates $`(u,v,\varphi )`$ as

$$x=\frac{R}{2}\mathrm{sin}u\mathrm{sinh}v\mathrm{cos}\varphi ,y=\frac{R}{2}\mathrm{sin}u\mathrm{sinh}v\mathrm{sin}\varphi ,z=\frac{R}{2}\mathrm{cos}u\mathrm{cosh}v.$$
(23)
and introduce the new fictive time scale $`\tau `$ satisfying

$$dt/d\tau =G(u,v)=r_Ar_B=R^2(\mathrm{cosh}^2v\mathrm{cos}^2u)/4.$$
(24)
Similar to the previous subsection, regarding the motion time $`t`$ and the negative total energy $`E`$ as generalized coordinate and generalized momentum respectively, then we can write the Hamiltonian function in extended space as

$$H^{}=HE0=H^{}(x,y,z,t;p_x,p_y,p_z,E)$$
(25)
Introducing the third category generating function $`F_3(p_x,p_y,p_z,E;u,v,\varphi ,t)`$, which satisfies

$$F_3=\frac{R}{2}\mathrm{sin}u\mathrm{sinh}v\mathrm{cos}\varphi p_x\frac{R}{2}\mathrm{sin}u\mathrm{sinh}v\mathrm{sin}\varphi p_y\frac{R}{2}\mathrm{cos}u\mathrm{cosh}vp_z+Et.$$
(26)
Thus momenta $`(p_u,p_v,p_\varphi ,E)`$ are related with old momenta $`(p_x,p_y,p_z,E)`$ by

$`p_u`$ $`=`$ $`{\displaystyle \frac{F_3}{u}}={\displaystyle \frac{R}{2}}\mathrm{cos}u\mathrm{sinh}v\mathrm{cos}\varphi p_x+{\displaystyle \frac{R}{2}}\mathrm{cos}u\mathrm{sinh}v\mathrm{sin}\varphi p_y{\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{sinh}vp_z,`$ (27)
$`p_v`$ $`=`$ $`{\displaystyle \frac{F_3}{v}}={\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{cosh}v\mathrm{cos}\varphi p_x+{\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{cosh}v\mathrm{sin}\varphi p_y{\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{cosh}vp_z,`$ (28)
$`p_\varphi `$ $`=`$ $`{\displaystyle \frac{F_3}{\varphi }}={\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{sinh}v\mathrm{sin}\varphi p_x+{\displaystyle \frac{R}{2}}\mathrm{sin}u\mathrm{sinh}v\mathrm{cos}\varphi p_y,`$ (29)
$`E`$ $`=`$ $`{\displaystyle \frac{F_3}{t}}=E.`$ (30)
And the regularized Hamiltonian is

$`K`$ $`=`$ $`{\displaystyle \frac{dt}{d\tau }}H^{}=G(u,v)(HE)0,`$ (31)
$`=`$ $`{\displaystyle \frac{1}{2}}[p_u^2+p_v^2+(\mathrm{csc}^2u+\mathrm{csc}h^2v)p_\varphi ^2]R\mathrm{cosh}v{\displaystyle \frac{(EH_i)}{4}}R^2(\mathrm{cosh}^2v\mathrm{cos}^2u).`$ (32)
where, $`H_i=z\epsilon (t)=(R/2)\mathrm{cos}u\mathrm{cosh}v\epsilon (t)`$. Obviously, $`p_\varphi `$ is a constant, which corresponds to the component of the angular momentum which along the $`zaxis`$. If $`p_\varphi `$ is equal to zero, the electron is constrained in a plane which can be chosen as $`y=0`$, corresponding to a two-dimensional motion. In the field-free case, $`\epsilon (t)=0`$, this two-dimensional model corresponds to the classical hydrogen molecular ion in ground-state. Equations of the motion for the ground-state hydrogen molecular ion interacting with laser field can be easily obtained from the regularized Hamiltonian, i.e.,

$`{\displaystyle \frac{du}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{K}{p_u}}=p_u,`$ (33)
$`{\displaystyle \frac{dp_u}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{K}{u}}={\displaystyle \frac{ER^2}{4}}\mathrm{sin}(2u)+\epsilon (t)(z{\displaystyle \frac{G}{u}}+G{\displaystyle \frac{z}{u}}),`$ (34)
$`{\displaystyle \frac{dv}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{K}{p_v}}=p_v,`$ (35)
$`{\displaystyle \frac{dp_v}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{K}{v}}=R\mathrm{sinh}v+{\displaystyle \frac{ER^2}{4}}\mathrm{sinh}(2v)+\epsilon (t)(z{\displaystyle \frac{G}{v}}+G{\displaystyle \frac{z}{v}}),`$ (36)
and

$$\frac{dE}{d\tau }=\frac{K}{t}=zG\frac{d\epsilon (t)}{dt},\text{ }\frac{dt}{d\tau }=\frac{K}{E}=G(u,v).$$
(37)
In the field-free case, defining

$`K_u`$ $`=`$ $`{\displaystyle \frac{1}{2}}(p_u^2+p_\varphi ^2\mathrm{csc}^2u)+{\displaystyle \frac{ER^2}{4}}\mathrm{cos}^2u,`$ (38)
$`K_v`$ $`=`$ $`{\displaystyle \frac{1}{2}}(p_v^2+p_\varphi ^2\mathrm{csc}h^2v)R\mathrm{cosh}v{\displaystyle \frac{ER^2}{4}}\mathrm{cosh}^2v.`$ (39)
thus $`K_u(=K_v)`$ is a constant of the field-free motion, they are related to $`\gamma `$ and $`\mathrm{\Omega }`$ by

$$K_u=K_v=\gamma R^2/4=ER^2/4+\mathrm{\Omega }/2.$$
(40)
Constants $`\mathrm{\Omega }`$ and $`\gamma `$ are first introduced by Erickson  and Strand  respectively, which have the following forms

$`\gamma `$ $`=`$ $`E2\mathrm{\Omega }/(mR^2),`$ (41)
$`\mathrm{\Omega }`$ $`=`$ $`\stackrel{}{L_A}\stackrel{}{L_B}+emR^2(\mathrm{cos}\theta _A\mathrm{cos}\theta _B).`$ (42)
Here, $`\stackrel{}{L_A}`$ and $`\stackrel{}{L_B}`$ are the angular momentum vector of the motion around nucleus $`A`$ and $`B`$ respectively, $`\theta _A`$ and $`\theta _B`$ are the angle from the vector $`\stackrel{}{r_A}`$ and $`\stackrel{}{r_B}`$ to positive $`zaxis`$ respectively, $`m`$ and $`e`$ are the mass and the charge of the electron respectively.

### C Action-angle variables and initial distributions

The action-angle variables for separated systems are defined as

$$I_i=\frac{1}{2\pi }p_i𝑑q_i,\text{ }\theta _i=\frac{}{I_i}p_i𝑑q_i=\frac{p_i}{H}\frac{H}{I_i}𝑑q_i.$$
(43)
The motion of the classical field-free hydrogen atom is periodic, then it need only a pair of conjugated action-angle variables

$$H_0=E_0=\frac{1}{2I^2}.$$
(44)
The corresponding angle is given by the Kepler’s equation

$$\theta =ue\mathrm{sin}u.$$
(45)
This means $`\theta `$ is the mean anomaly of the free orbits. The eccentric anomaly $`u`$ can be obtained from $`r=a(1e\mathrm{cos}u)`$. Here, $`r`$ is the distance from the electron to the origin, $`a`$ is the instantaneous semimajor axis which satisfies $`a=1/(2E_0)`$, $`e`$ is the eccentricity of the free orbits satisfying $`e=\sqrt{2E_0L^2+1}`$, $`L`$ is the total angular momentum.

The motion of the ground-state hydrogen molecular ion $`H_2^+`$ is quasi-periodic, it need two pairs of action-angle variables. With elliptic integrals, action $`I_u`$ can be expressed as follows.

$$I_u=\{\begin{array}{cc}\frac{\sqrt{8K_u2E_0R^2}}{\pi }F_1(\frac{\pi }{2},\sqrt{1+\frac{4K_u}{E_0R^24K_u}}),\hfill & for\text{ }K_u>0,\hfill \\ \sqrt{E_0R^2/\pi },\hfill & for\text{ }K_u=0,\hfill \\ \frac{\sqrt{2E_0R^2}}{\pi }[F_1(\frac{\pi }{2},\sqrt{1\frac{4K_u}{E_0R^2}})\frac{4K_u}{E_0R^2}F_2(\frac{\pi }{2},\sqrt{1\frac{4K_u}{E_0R^2}})],\hfill & for\text{ }K_u<0.\hfill \end{array}$$
(46)
where, the first category elliptic integral $`F_1(\phi ,k)`$ and the second category elliptic integral $`F_2(\phi ,k)`$ are in the form of

$`F_1(\phi ,k)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\phi }{}}}\sqrt{1k^2\mathrm{sin}^2x}𝑑x,`$
$`F_2(\phi ,k)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\phi }{}}}1/\sqrt{1k^2\mathrm{sin}^2x}𝑑x.`$
Generally, a single trajectory lacks the inversion symmetry of the real physical systems. So the spectrum obtained from a single trajectory exhibits unphysical even harmonics. A natural way to remedy the unphysical even harmonics is to consider an ensemble of trajectories, evolving from an initial microcanonical distribution with inversion symmetry. For a chaotic system, the initial distribution can be generated with Monte-Carlo method. However, for an integrable system, the distribution generated by Monte-Carlo method does not possess of ergodicity. We find that the points on the same equienergy surface generated by action-angle variables with regular steps possess of good ergodicity. To reconstruct the inversion symmetry, there must exist pairs of $`(x_0,y_0,z_0;p_{x_0},p_{y_0},p_{z_0})`$ and $`(x_0,y_0,z_0;p_{x_0},p_{y_0},p_{z_0})`$ in the initial distribution. For the regularized model for hydrogen atom, it corresponds to pairs of $`(u_0,v_0,\varphi _0;p_{u_0},p_{v_0},p_{\varphi _0})`$ and $`(u_0,v_0,\varphi _0+\pi ;p_{u_0},p_{v_0},p_{\varphi _0})`$. And for the regularized hydrogen molecular ion $`H_2^+`$, it corresponds to pairs of $`(u_0,v_0,\varphi _0;p_{u_0},p_{v_0},p_{\varphi _0})`$ and $`(u_0+\pi ,v_0,\varphi _0;p_{u_0},p_{v_0},p_{\varphi _0})`$.

## III Numerical simulation

The numerical computational procedure is based upon the classical trajectory Monte Carlo (CTMC) method \[20-21\]. CTMC simulation procedure involves three stages, (i) choice of initial conditions, (ii) numerical integration of equation of motion, and (iii) categorization of each trajectory as excitation, charge transfer or ionization. In the process of numerical integrating, numerical accuracy and computing time are two primary aspects that must be considered. We use the fourth-order Runge-Kutta method with variable steps to perform the numerical calculation. Note also that computer can not deal with singularity, which corresponds to electron-nucleus collision, to overcome this difficulty, we have implemented regularization. With complete regularization, numerical simulation can be established with required precision before, at, and, after collision successfully.

To obtain the harmonic spectra, our procedure also calculate the averaged dipole moment of the excited trajectories with pairs of inversion symmetric initial conditions in the same distribution. Having determined the actual trajectories of the electron, one can easily obtain the component of the averaged dipole moment, which along the laser polarization direction. Then the harmonic spectra of the driven dipole is straightforwardly obtained from it’s power spectra

$$D(\omega )=\underset{t+\mathrm{}}{lim}\frac{1}{t}\left|\underset{0}{\overset{t}{}}\mu (\tau )e^{i\omega \tau }𝑑\tau \right|^2.$$
(47)
where, $`\mu (\tau )`$ is the component of the averaged dipole moment along the laser polarization direction. In our calculation, the electric fields of the ultrashort laser pulses are chosen as

$$\epsilon (t)=\{\begin{array}{cc}E_M\mathrm{sin}^2(\omega _Lt/40)\mathrm{sin}(\omega _Lt),\hfill & for\text{ }0t40\pi /\omega _L,\hfill \\ 0,\hfill & otherwise.\hfill \end{array}$$
(48)
The maximum electric field strength $`E_M`$ is related to the laser intensity $`I=\sqrt{\epsilon _0/\mu _0}E_M^2/2`$, and $`\omega _L`$ is the angular frequency, and period $`T_L`$ is equal to $`2\pi /\omega _L`$.

The time-dependent dipole moment $`\mu _z(t)`$ of a single trajectory along the laser polarization direction sensitively depends on the laser intensity. Fig.1 shows the time evolution of the dipole of the hydrogen molecular ion with initial energy $`E_0=1.1034`$ $`hatree`$, internuclear distance $`R=2.00`$ $`bohr`$, laser wavelength $`\lambda =600`$ $`nm`$ and different laser intensity. And some dipoles of the hydrogen atom and their power spectra are presented in Fig.2, with initial energy $`E_0=0.5`$ $`hatree`$, laser wavelength $`\lambda =532`$ $`nm`$ and different laser intensity. With the increasing of the laser intensity, the oscillations of the dipole moments are modulated gradually, and it follows a regular pattern with both high and low frequency components. Such patterns are independent of initial conditions. When the laser intensity is large enough, ionization takes place, which strongly modifies the subsequent time evolution of the dipole moment.

The power spectra obtained from a single trajectory of the laser driven hydrogen molecular ion, which with initial energy $`E_0=1.1034`$ $`hatree`$, internuclear distance $`R=2.00`$ $`bohr`$, laser wavelength $`\lambda =600`$ $`nm`$ and different laser intensity, are presented in Fig.3. The first one corresponds to the field-free case. For the hydrogen atom, it possesses of only regularly decreasing peaks locating at multiples (harmonics) of Kepler’s frequency, that is, $`\omega _n=n\omega _0`$, for $`n=1,2,3\mathrm{}`$, which manifests that the free motion is periodic. While for the hydrogen molecular ion, the free motion is quasi-periodic, it appears regularly decreasing peaks locating at two different characteristic frequencies and their combinations, i.e., at $`n\omega _{01}+m\omega _{02}`$, for $`n=0,1,2,3\mathrm{}`$, $`m=0,1,2,3\mathrm{}`$ and $`n+m>0`$. For the laser driven systems, the peaks of harmonic spectra depend strongly on the laser intensity. In addition to the origin peaks, a dominant line at the laser frequency $`\omega _L`$ appears, which is the Rayleigh component in the light scattered by the atoms and molecular ions. High-order harmonics, which consist of both odd and even components, do appear in the spectra, and their orders and strengths increase with the increasing of the laser intensity. For a low laser intensity, i.e., in the perturbative regime, the characteristic peaks dominate the power spectrum. For a proper laser intensity, the plateau structure containing both odd and even harmonics appears. For a larger laser intensity, because of the presence of chaos, the spectra become noisier even the motion is bound. For a strong enough laser intensity, the motion becomes unbound and the corresponding spectrum is dominated by a very noise background, which is generated by the ionized electrons, the only line surviving being the one at the laser frequency, which corresponds to the light scattered by the asymptotically free electron Thomson scattering.

To eliminate the unphysical even harmonics, averaging an ensemble of trajectories evolving from an inversion symmetric distribution is necessary. The time evolution of the averaged dipole, which evolve from a microcanonical ensemble of 5000 trajectories, are presented in Fig.4, the first one corresponds to the hydrogen atom with initial energy $`E_0=0.5`$ $`hatree`$, laser wavelength $`\lambda =532`$ $`nm`$ and laser intensity $`I=5.0\times 10^{14}`$ $`W/cm^2`$, the second one corresponds to the hydrogen molecular ion with initial energy $`E_0=1.1034`$ $`hatree`$, internuclear distance $`R=2.00`$ $`bohr`$, laser wavelength $`\lambda =600`$ $`nm`$ and laser intensity $`I=1.0\times 10^{14}`$ $`W/cm^2`$. Comparing with the evolution of the dipole of a single trajectory, one can easily find that the oscillation of the averaged dipole is smooth and it globally follows the laser oscillation.

The spectra obtained from an ensemble of trajectories are showed in Fig.5 and Fig.6. Fig.5 is the harmonic spectra of the ground-state hydrogen atoms interacting with the laser pulses with $`\lambda =532`$ $`nm`$ and different laser intensity. Fig.6 is harmonic spectra of the hydrogen molecular ion with initial energy $`E_0=1.1034`$ $`hatree`$, internuclear distance $`R=2.00`$ $`bohr`$, laser wavelength $`\lambda =600`$ $`nm`$ and different laser intensity. As a consequence of averaging process, the unphysical even harmonics are remedied really. Within a proper laser intensity range, the plateau structure that only possesses of odd harmonics appears. As pointed out previously, at a higher intensity, the spectra become noisier even the ionization does not happen because of the effects of chaos. This indicates that the chaos cause the noise of the harmonic spectra. When the laser intensity is high enough, the ionization takes place, thus the noise background conceals the high-order harmonics. This means that the onset of ionization follows the onset of chaos and the ionization suppresses the harmonic generation.

## IV Summary and discussion

In summary, within the Born-Oppenheimer approximation and using the classical trajectory method,we have calculated the high-order harmonic generation spectra of the hydrogen atom and the hydrogen molecular ion interacting with ultrashort intense laser pulses. The other multi-photon phenomena, such as multi-photon ionization and above threshold dissociation, can also be simulated. In our subsequent calculations, the dynamics of the electron is investigated by numerical integrating the equations of motion using regularized coordinates. To eliminate the unphysical even harmonics of a single trajectory, averaging over an ensemble of trajectories evolving from an initial microcanonical distribution with inversion symmetry is necessary. Such distribution is constructed using action-angle variables. A plateau structure in the spectra with only odd harmonics is observed within a proper laser intensity range of about $`10^{14}`$ $`W/cm^2`$. From our numerical results, we observe that the high order harmonics are cut off at a special order harmonic. At higher laser intensities, chaos introduces noise into the spectra even though the motion is still bound. Finally as the intensity is further increased, ionization takes place, and the harmonics disappear.

These results are qualitatively consistent with recent quantum calculations  and experimental observations \[1-5\], but the cutoff order $`N_m`$ of the plateau structure is not precisely consistent with formula $`N_m=(I_p+3.17U_p)/\mathrm{}\omega __L`$, here, $`I_p`$ is the ionization potential and $`U_p=e^2E_M^2/4m_e\omega __L^2`$ denotes the quiver energy or the ponderomotive energy of an electron. As an example, when the laser intensity $`I=10^{14}`$ $`W/cm^2`$ and wavelength $`\lambda =600`$ $`nm`$, the ionization potential $`I_p`$ of the ground-state hydrogen molecular ion is $`1.1034`$ $`hatree`$ $`(29.77eV=14.50`$ $`\mathrm{}\omega __L)`$, the quiver energy $`U_p=3.36eV`$ $`(=1.63`$ $`\mathrm{}\omega __L)`$, then $`N_m=19`$, and when $`I=7.5\times 10^{13}`$ $`W/cm^2`$, $`N_m=18`$; however, from our simulation, the harmonic plateau are both cut off at $`17`$, and for very large intensity (above $`10^{15}`$ $`W/cm^2`$) the high-order harmonics are concealed by the noises. To obtain quantitative results, we have to integrate the time-dependent Schr$`\stackrel{..}{o}`$dinger equation, this can be realized with the split-operator method . For the hydrogen molecular ion, within BOA, a proper internuclear distance $`R`$ will enhance the high-order harmonic generation, and it will be interesting to go beyond the Born-Oppenheimer approximation to investigate what further interesting insights can be obtained when the nuclear motion is taken into account \[23-24\].

In our model, we only consider the non-relativistic case with dipole approximation. When the laser is sufficiently intense, photoelectrons of relativistic energies can be produced, necessitating a fully relativistic treatment . The dipole approximation is no longer valid, and the magnetic field is not only present, but acquires an importance similar to that of the electric field. And even before the appearance of the relativistic photoelectrons, the effects of the magnetic field may be very important too . Due to the wiggly motion and the acceleration of the electron near the outermost turning points induced by the magnetic field, the cut-off order of the harmonic plateau maybe higher.

Acknowledgment

The work is supported by the National Natural Science Foundations of China under Grant No. 19874019, 19904013 and 1990414. The author Chaohong Lee thanks very much for the help of Dr. Haoseng Zeng, Dr. Ming He and Dr. Zongxiu Nie.


text:
# Electronic structure and magnetism in slightly doped SrB6.

## I Introduction.

The structure of hexaborides, XB6, is simple cubic, with one large atom (X = Eu, Ca, Sr or Ba) in the corner, and a B6 octahedron in the center. It can be viewed as the CsCl-structure, where the B6-octahedron would be a pseudo-atom in the center and the X-atom in the corner is the second atom. The internal parameter, $`z`$, is a measure of the size of the octahedron. It is the shortest distance from one of the B-atoms at the edge of an octahedron to the cube surface. Its value is about 0.2 of the lattice constant, but it varies from system to system. The hexaborides have been studied mostly because their intriguing properties showing signatures of both metallic and semi-conducting behavior . Recent experimental works on La-doped hexaborides like LaxCa(1-x)B6 show unexpectedly a very weak ferromagnetic state that persists up to large temperatures, and for very low $`x`$, about 0.005 . Independent measurements confirm these findings in part, but they also show that the results can span from weak ferromagnetism to strong paramagnetism and that Al impurities may conceal some of the intrinsic properties .

Ferromagnetism in the La-doped hexaborides XB6, X=Ca,Sr or Ba, is surprising because the density-of-states (DOS) at the Fermi energy, $`N(E_F)`$, and the Stoner factor, $`\overline{S}`$, are expected to be small in this family of compounds. Band calculations  in the local density approximation (LDA)  on SrB6 and CaB6 confirm that the DOS is low, with $`E_F`$ falling in a valley of the DOS, and spin-polarized calculations for EuB6 show a very weak spill-over of the spin density from the localized Eu-f spin to the rest of the valence band. The exact behavior of the band overlap at the X-point has been monitored in several band calculations . It has been found that this overlap, which determines if the undoped system is metallic or semi-conducting, depends on the structural parameter $`z`$. In fact, the overlap is more sensitive to $`z`$-variations than to variations of the lattice constant . A recent band calculation  determined the fine details of the bands near $`E_F`$ and the Fermi surface (FS), and concluded that they are in good agreement with de-Haas-van-Alphen measurements . An experimental determination of the band dispersion the near the X-point, by photoemisson , is in fair agreement with the theoretical predictions for the existence of a X-pocket in the bulk, but a different interpretation is made for bands in the surface region.

Electronic structure calculations based on the local (spin) density approximation, either by spin-polarized calculations or by calculations of the Stoner factor $`\overline{S}`$, are in general able to detect if a system is ferromagnetic or not . Because the Stoner factor for LaxSr(1-x)B6 is far below the limit for magnetism , it is tempting to suggest that the observed magnetism is not a band magnetism, but something else. Electron-hole interaction for an asymmetric position of EF relative to the gap has been proposed to lead to excitonic magnetism . Spontaneous polarization of a dilute electron-gas has also been suggested . However, for Ce it is found that the Stoner factor underestimates the tendency towards magnetism, because the stabilization energy for a magnetic state comes also from the Coulomb energy in addition to the exchange energy . The effects can be larger in a compound as long as the different sub-bands and the ionicity of the system show favorable conditions. In a recent work we showed that LaSr7B48 is likely to be such a system . Here we extend these calculations to larger unit cells, to different types of dopants, and give a more complete discussion of the results.

## II Theory.

### A Band method.

Electron energies are calculated by the linear muffin-tin orbital (LMTO) method  using the LDA. The internal parameter $`z`$, is in these calculations 0.215 of the lattice constant. The lattice constant is fixed at 7.93 a.u. in each calculation. Additional ’empty’ spheres are inserted at the cube edges, between the Sr-atoms, in order to account for non-sphericity of the potential in the open part of the structure. Calculated results are presented for two unit cells: One supercell of 8 elementary cubes (2 cube lengths in x,y and z-direction), containing totally 80 sites, 56 atoms and 24 empty spheres. A second set of calculations was made using a larger super cell containing 27 elementary cells (3 cube lengths in each direction), with totally 270 sites (189 atoms plus 81 empty spheres). These two cells allow to study doping concentrations $`x`$ of 0.125 and 0.037, respectively. Nonequivalent sites are treated independently. One basis set for the 80-site cell unit cell includes s,p,d for the atoms, and s,p for empty spheres. The results from these calculations can be found in ref. . A reduced basis set for the large 270-site calculations includes s and p states for the B-atoms and only s-states for the empty sites in the direct basis. We include the B-d and empty-p ”tail” states in the reduced basis calculations. This means that the d-states for B atoms, for example, are constructed from tail states coming from neighboring sites. This gives a good description of tail states like B-d and p-states on the empty sites, close to what was obtained from the large basis, without having to increase the size of the eigenvalue problem. Intrinsic ”on-site” states, like the high energy La-f state or B-d states, can not be constructed from tail states in this way, but these states are unoccupied and not important in this study. Test calculations for the smaller of the cells using this smaller basis showed very similar results as the ones from using the large basis, when the linearization energies (E) of the La-band are increased upwards about 0.2 Ry.

In the 270-site calculations, where we employed the same basis (with tail states) and the same E’s, which in the small basis 80-site calculation reproduced well the results of the older large basis calculation. In this way we believe that the results for the large 270-site cell are realistic despite the use of a minimal basis set. No lattice relaxation around the impurities is considered. The number of k-points varies from 10 to 20 irreducible k-points for the unit cell of 80 sites, and from 4 to 10 for the 270 site cell.

### B Stoner models.

The original Stoner model  can be applied together with ab-initio paramagnetic band results in order to study the spin susceptibility or the onset of magnetism . An applied exchange splitting, $`\xi `$, leads to a transfer from minority spins to majority spins, and a loss in kinetic energy $`\mathrm{\Delta }K=N\xi ^2`$, where $`N`$, the DOS at $`E_F`$, is assumed to be constant. The same spin transfer leads to a gain in exchange energy $`\mathrm{\Delta }E_x=N^2I_s\xi ^2`$, where $`I_s`$ is the exchange integral . If a partial DOS is not constant near $`E_F`$, there will be a possibility of charge transfers, $`\mathrm{\Delta }q_\xi `$ and $`\mathrm{\Delta }Q`$, as function of $`\xi `$ and temperature, which make additional contributions to the total energy. It is possible to identify one sub-band, the La-impurity (mainly d) band, which has a large derivative of its DOS, $`N_{La}^{}(E_F)`$. The total DOS at $`E_F`$ is $`N=N_{La}+N_v`$, where $`N_v`$ is the rest of the DOS. By neglecting all other derivatives one can make a model of the charge transfers as function of $`\xi `$ and $`k_BT`$. By taking into account the energy associated with the charge transfer, $`U_0`$, one obtains a modified Stoner criterion. This is made in detail previously  and the new Stoner criterion is:

$$\overline{S}_U=NI_s+U_0N_{La}^{}N_v/N^21$$
(1)
This modified Stoner criterion is useful for an understanding of why a system can be pushed towards magnetism by Coulomb energies. But, as will be discussed in section 3.2, the effects are too strong to use eq. 1 for a quantitative prediction of a magnetic state.

## III Results.

### A Bandstructures.

The bandstructure for undoped SrB6 agree well with the one in ref. . The Fermi energy falls in a low DOS region between the B-p and Sr-d band regions, where a small gap can appear. The B octahedron behaves partly as one identity, as a large pseudo atom. This is seen in a narrow band about 0.9 Ry below EF. There is only one such band per basic unit cell, although in reality it is of B s- and p-character. This band is separate from the other B s- and p-bands. Its presence could be viewed as one single s-band from each B6 pseudo atom within a CsCl-structure, where Sr is one atom and B6 the other one. However, this is far below EF and this way of looking at the B-octahedron is probably not useful for the bands near EF or for the physical properties of doped hexaborides.

The effect of doping, where one Sr is exchanged with another atom, show some unexpected result depending on the dopant. Fig 1 shows the DOS for undoped, for La-doped and In-doped SrB6, all for $`x`$=0.125 (small 80-site cell, small basis). La has a d-band coming in just above EF. In this energy range there are very few B and Sr bands for the La-d band to hybridize with. The La-d band can hybridize between La neighbors if the La concentration is large, but if there are only few La atoms, the La-d band will localize and becomes narrow. Indium has one additional valence electron compared to Sr, exactly as La has. But the additional electron in In is a p-electron, and in contrast to the d-band for La, the In-p band hybridizes well with the B-bands below EF. The new In-band below EF can accommodate two electrons, but In has only one additional electron compared to Sr. The result is that EF will move downwards compared to the rest of the bands. Thus, as seen in fig. 1, the effects of In and La doping are opposite. Furthermore, with La there is a large increase of the local La-DOS, while the DOS on the In-site and a Sr-site are very similar. The resulting DOS from In doping is very similar to the one from Al (or even B) doping, at least as long as lattice relaxation can be ignored. It has been reported that Al can penetrate the hexaboride material when the aluminum flux-growth technique is used to prepare single crystals . The calculated results for single In or Al impurities shown here, indicate that the effects of Al contamination should be very different from those of La impurities. The effects could even appear quite undetected compared to undoped hexaboride, since the main result is a rigid-band like downward shift of EF on a uniform DOS. Theories of excitonic magnetism are based on rigid-band assumptions of the undoped band structure . Such theories should therefore apply to Al (or In) doping, but not easily for La-doping, since the latter is not rigid-band like. However, according to our results, the effect of the p-electron doping in this material is oddly transformed into hole doping. The charge distribution is such that about one half of the additional electron from In stays on the site, a bit more than one half is added close to the impurity. At the most distant sites there is slightly less electronic charge compared to the values for the undoped case, which presumably is an effect of the shift of EF.

The paramagnetic (non-spinpolarized) DOS from the calculation for LaSr7B48 is shown in fig. 1. The La band forms a separate band near $`E_F`$ below the rest of the Sr-d band, making the real space distribution to be localized near the La impurity from this band. The total LaSr7B48 DOS at $`E_F,N(E_F)`$, about 48 states/cell/Ry, is concentrated on the single La site, with $`N_{La}(E_F)`$ 19 states/cell/Ry, cf Table 1. On the remaining 55 sites, each Sr has on average 1.0 and each B roughly 0.4 states/cell/Ry, respectively. By comparison, a rigid-band model obtained by adding one electron to the band structure without the La impurity gives very different values: Of the total DOS, 32 states/cell/Ry, each Sr has 1.4 and each B 0.25 states/Ry. Clearly, La behaves like a localized impurity site with a very large local $`N_{La}(E_F)`$, much larger than on Sr. The DOS values are included in Table 1.

The additional valence electron remains located on the La. By comparing the charges in the undoped and doped supercell one finds that the remaining Sr$`{}_{7}{}^{}B_{48}^{}`$-spheres receive no additional charge, cf. Table 2. The difference in charge between the two systems is never larger than 0.01 electron per Sr or B site. The derivative of the DOS becomes very large, about 1400 states/cell/Ry2, from which 700 comes from the La site, compared to about 30 for each Sr site when the rigid-band model is applied to the undoped bandstructure. This asymmetric and uneven distribution of the DOS near $`E_F`$ is a condition for real space charge transfers as T increases and as an exchange splitting is imposed. The transfer goes in both cases to the La site from the rest of the system.

Next, for the calculations with the large 270-site cell only the case with a La-impurity was considered, LaSr26B162. The total and partial La-DOS near EF are shown in fig 2, together with the corresponding DOS from the 80-site calculation. Both results are calculated with the same basis. The DOS functions from the two calculations are divided by 27 and 8, respectively. This makes the units of the DOS states per formula unit (f.u. of SrB6), and the two cases can be compared on the same diagram. The number of k-points, 20 for the small and 10 for the large cell, reflects a larger density of points in the reciprocal space for the 270-site calculation. It is seen that the La-band is more separated from other bands and more localized in the large cell. The value of the total N at EF (6 states per Ry and f.u.) is quite similar in the two cases, but the local DOS on the La-site is much larger in the largest cell. The 270-site cell has a large number of sites far from the impurity with low DOS, almost as for the undoped material. The number of such sites is comparatively larger than for the 80-site cell. But the very large DOS at the impurity site compensates for the low DOS elsewhere, and the average DOS per f.u. becomes similar for the two cells. The local La-DOS, which in the small cell is about 19, is 55 in the large cell, cf Table 3.

A comparison of the DOS near EF with the ones from the small cell, reveals more rapid variations of the DOS structures. The La-band is narrower, with the total DOS going from zero to very large values within a much narrower energy interval than in the case of the 80-site cell. The effect is larger DOS derivatives, mainly steming from La, but also larger sensitivity of the charge transfers and their T-dependencies. The amount of charge transfer as function of $`k_BT`$ (and $`\xi `$) is typically 3-5 times larger within the large cell, but the non-linearities are large and a small shift of EF can make a large difference. It is therefore not meaningful to quote values of $`N^{}`$ in this case. At the most distant (from the La) B sites, the local B-DOS is small, but still about 50 percent larger than what is the case for the undoped DOS, when a rigid band shift has been done to account for one additional electron. For the distant Sr-sites, the local DOS tend to be close to the value of the undoped DOS with the proper rigid-band shift of EF, about 0.7 per Ry. The first 1-2 layers of Sr and B surrounding the La impurity, have approximately doubled their local DOS compared to the undoped values.

Despite the large number of atoms separating the La in the largest cell (188 atoms, 3 types of Sr and 9 types of B), and despite the fact that the additional La electron stays localized near the impurity, there is still an influence on the electronic structure far from the impurity. Their DOS values are larger than the ones of the undoped DOS, which are obtained after a rigid-band shift of EF to account for an additional electron. An even lower impurity concentration seems to be required for finding the electronic structure far from the impurity, almost ”insulating” like for the undoped case. In such a case, the interaction between local moments would effectively be cut-off.

### B Magnetism.

We first present the results using the 80-site cells, mostly using the small basis set as explained above. Spin-polarized calculations for undoped SrB6 do not give stable moments of significant size. An applied field is used to get a moment of about 0.5 $`\mu _B`$ per cell, but this moment decreases steadily during the iterations when the field is removed. The same situation holds for InSr7B48. A number of calculations were also attempted where one Sr was replaced with an ’empty’ site. This was motivated by the reports that a kind of ”self-doping”, by Sr vacancies might be present because some seemingly undoped samples show the weak magnetism . The unpolarized calculations for the cases with a Sr-vacancy showed the effect of charge transfer as function of T, so that as T increased there was a larger charge on the remaining Sr sites. However, the moments always decreased (to below 0.005 $`\mu _B`$) in the spin-polarized calculations, also for the case where a slight B-relaxation towards the missing Sr-site was attempted. Thus, in the following we discuss only La-doping, since undoped, In- (or Al-) doping, and Sr-vacancy ”self-doping” all failed to show magnetism.

The 80-site cell calculation of LaSr7B48 using 20 k-points and the large basis, finds a small moment of 0.10 $`\mu _B`$/cell, which is about 10 mRy/supercell lower in total energy than the paramagnetic solution . If fewer k-points are used, one finds that the (paramagnetic) DOS and its energy derivative near $`E_F`$ are larger, cf Table 1. One case with the small basis is also shown. The moments are not quite as large as with the large basis, but qualitatively the results are the same.

A number of different calculations with the 80-site cell and the small basis are done where the linearization energies are varied in order to obtain different values for N(EF) and its derivative. All these different cases are not presented in detail, but the trend is that when the local La-DOS and its derivative at EF are large, then the moments will be large (larger than about 0.1 $`\mu _B`$ per cell) in the spin-polarized calculations. Thermal broadening will usually oppose a large moment. But there are a few exceptions, when the energy derivative of the DOS is varying around EF, showing that the sharpness of the DOS is not always close to EF. Temperature dependent charge transfers also show deviations from a simple T2 behavior. The moment is always largest on La, with more than 50 percent of the total moment on and near the La impurity. The moment on B-sites far from the impurity is negligible. This fact and the uneven distribution of local DOS show that a rigid-band picture of dilute La-doped hexaborides does not apply. The bands near $`E_F`$ are modified by the spatially non-uniform perturbation of the potential around the impurity site.

The calculated Stoner product $`NI_s`$ is 0.21 , $`N_{La}^{}`$ 730 per Ry$`{}_{}{}^{2}`$ cell, N $`48`$ and $`N_v29`$ (per Ry $``$ cell), to give $`\mathrm{\Delta }q_\xi 440\xi ^2`$ and $`\mathrm{\Delta }Q=720(k_BT)^2`$, with $`\xi `$ and $`k_BT`$ in Ry. With $`U_0`$ calculated from T-dependent variations of the total energies to be about 1.5 Ry/el./La, the Coulomb term in eq. 1 becomes very large, about 14. However, even if the Coulomb energy represents an important correction in this case, there are several facts indicating that other energies are involved, and that the response of the system (during a self-consistent, spin-polarized calculation) will prevent the system to absorb all of the energy. First, an exchange splitting $`\xi `$ which is not uniform over all sites will modify the energy. The spin-polarized calculation shows that $`\xi `$ is finite (3-5 mRy) only close to La. Second, although the amount of charge transfer due to $`\xi `$ and $`k_BT`$ is similar, there are long tails of the Fermi-Dirac distribution which make the radial dependencies of $`\mathrm{\Delta }q_\xi `$ and $`\mathrm{\Delta }Q(T)`$ different and make $`U_0`$ non-transferable from one case to the other. Third, the effective derivative of the La-DOS, $`N_{La}^{}(ϵ,T)`$, (taken as $`(\frac{df}{dϵ})N_f^{}(ϵ)𝑑ϵ`$), is T-dependent, because the bare DOS does not have the same large derivative over a wide energy range. The mechanism of charge transfer and magnetism is therefore expected to depend on T. This is confirmed in the spinpolarized results. Usually, by raising $`k_BT`$, the magnetic moment goes down. But a few cases are found where the moment is largest for $`k_BT`$ about 3 mRy. The reason is that the sharpest DOS structure is a bit beside the position of EF in that case.

Calculations of total energies converge slowly and a comparison between paramagnetic and spin-polarized results are not made in all cases. From the converged cases it is seen that the total energies for magnetic solutions are lower than the non-magnetic ones. This difference is largest when the moments are largest, 10-14 mRy/cell for a total moment of 0.10-0.19 $`\mu _B`$. For the cases with moments smaller than 0.05 $`\mu _B`$, the difference is only 0-2 mRy.

Two sets of spin-polarized calculations were done for the 270-site cells, both with 4 and 10 k-points. One used the same linearization energies, which for the 80-site cell reproduced well the result with the large basis. It is the paramagnetic DOS functions from these parametrizations which are shown in fig. 2. The DOS values at EF and magnetic results are displayed in Table 3. The other set use different linearization energies, which for the paramagnetic 80- and 270-site calculations lead to increased DOS and derivative of the DOS. The spin-polarized results using 10 k-points show a much larger moment for the latter case compared to the first set. The reason is that there is a region of a rather flat (or less steep) DOS close to EF in the first case, as can be seen in fig. 2. The bold line of the DOS from the 270-site cell is seen to follow the thin line of the 80-site cell, but only within a few mRy of EF. Further below EF the curve is again steep, as is the case for the other set of calculations for the 270-site cell. The result is that the charge transfer is lower in the former set of calculations. However, this difference between the two sets of calculations diminishes at larger T when the steep region of the DOS a bit aside of EF becomes involved. This is consistent with the T-dependence of the moment, since the largest moment in the first set of calculations shown in Table 3, is when $`k_BT`$ is 4 mRy. At low T ($`k_BT`$ 0.5 mRy), the moment decreases further. The T-dependence of the moment is in general opposite for the other cases, meaning that the moments become very small when $`k_BT`$ 3 mRy. For the other set and in most cases with the 80-site calculations the moments are largest for the smallest $`k_BT`$ 0.5 mRy.

When the virtual crystal approximation is applied to this case, to remove 0.5 electrons from the system (and 0.5 nuclear charges from the Sr sites), one finds that EF falls closer to the steep left-hand side of the La-DOS, cf. fig 2. The objective with this model calculation is to see if the new conditions make the moment to increase, as was expected from the previous discussion. It should also approach the situation with two-component doping, where Al-doping (and a rigid-band shift of EF) in addition to La-doping could influence the magnetic state. The result is that the spin-polarized results converge towards a larger moment, about 0.09 $`\mu _B`$, but now only for low T. At large $`k_BT3`$ mRy, the moment shows a steady decrease. These T-dependent results corroborate the finding that there is a close connection between sharp DOS structures, charge transfer and magnetic moments.

The results in Table 1 and 3 show that the moments per impurity are comparable for the small and the large cell, about 0.1 $`\mu _B`$ per cell. The correlation between large derivative of the DOS and the moment, is seen in the results, both from the variation of the number of k-points and from the different choice of linearization. The largest moment, 0.25 $`\mu _B`$/La is only slightly larger than the largest moment in the 80-site calculation, 0.19 $`\mu _B`$/La. This is despite the much larger derivative of the DOS from the results with 270 sites. As was discussed above, the reason could be that the large cell has a comparatively larger undoped bulk, which opposes polarization.

### C Discussion.

The results in ref.  for an 80-site cell, showed a stable ferromagnetic state, and the moment remains centered around the La impurity. Together with the paramagnetic results, one can conclude that the magnetism is correlated with the derivative of the DOS on La, and that the magnetic state is partly stabilized by a charge localization on the La impurity. This correlation becomes more evident in the continued calculations using smaller basis, which are modified in order to create different conditions for large or small charge transfer effects. The results from the 270-site calculations show that the conditions for charge transfer to the La impurity are reinforced, because the La-d band becomes even more narrow and localized due to the increased La-La distance in the large cell. The magnetic moments per La-impurity are comparable in the two cells though, probably because of the larger ”bulk” between the impurities, which is undoped and counteracts magnetism better in the large cell. On the other hand there is no magnetic state for the In or Al impurity, for the Sr-vacancy or for the undoped SrB6. This is despite a certain pre-condition for a charge transfer from B towards remaining Sr in all these cases, but apparently it is too weak to produce an itinerant magnetic state, at least in the 80-site cell. The weak transfer towards Sr is present even in the La-doped cases, and therefore it reduces the (larger) transfer towards La. This competition could reduce the amplitude of the La-moment and even prevent magnetism if it became larger. One could speculate that Al substitutions on Sr sites play a role of a charge reservoir in this situation. Namely, as was noted earlier, an Al atom on a Sr-site makes no large effects on the DOS except for a rigid-band like downward shift of EF. No explicit calculations of having both La- and Al-substitutions were done, but taken separately these results show impurity band and rigid-band behavior, respectively. A varying Al content in a La-doped material should therefore modify the position of EF relative to the narrow La-d band, and create more or less optimal conditions for the charge transfer around the La-impurity. By tuning with Al-doping to find the largest DOS derivative on the La-band, one should be able to optimize the magnetic moment.

The non-rigid band like electronic structure around the dilute La impurity shows up in the distribution of the ferromagnetic moment. According to these results, it is striking that La and the nearest sites takes the major part of the total moment. B-sites remain essentially non-polarized. The La-site acts as an impurity and will be a perturbation for excitons. Ferromagnetism has been proposed to be caused by the excitonic mechanism , since the undoped band structure show some favorable conditions for this. However, it is seen that with La-doping there are large modifications of the undoped bandstructure, and a rigid-band assumption is not at all appropriate. Rigid-band like conditions occur for In-, Al-, to some extent for ”vacancy”-doping, with EF lower than for the undoped case. Theories of spontaneous ferromagnetism of the dilute electron gas have also been proposed . However, an unrealistically large electron-gas parameter $`r_s`$ is necessary for this ($``$ 30), whereas the calculated $`r_s`$ is smaller than about 1.5.

Some experiments can be suggested to verify if the La-band is localized and if the moment is concentrated close to the La. Nuclear magnetic resonance experiments can detect a magnetic moment. But the magnetic moment is small and for a total moment of 0.19 $`\mu _B`$ one finds a La contact field of less than about 40 kG. The spin-splitting is not exceeding 5 mRy within the La valence and core levels, which is small to be observed directly by spectroscopic methods. Furthermore, the low La concentration for the real magnetic system might be a problem for the intensity. The moment is mainly localized to the impurity, with weak interactions with the exterior. This situation resembles that of a Kondo impurity, but due to the small moment it is questionnable if this can produce the typical manifestations on specific heat and resistivity at low T. Instead of looking for the moment directly, one can propose to use spectroscopic methods like photoemission, soft X-ray emission spectroscopy (XES) or soft X-ray absorption (SXA), to determine the degree of localization of the La-band. This can be done at larger La concentrations, before magnetism sets in. Such measurements could reveal if the DOS shows a localized La-band similar to the DOS calculated here for $`x0.125`$ or even more at $`x0.037`$, or if the DOS is more like a rigid-band image of the undoped DOS. However, the low concentration of La can be a problem in some cases and lead to a inconclusive result. Soft x-ray photoemission (XPS), for instance, gives a often a good picture of the total DOS modulated by matrix elements, but if the La concentration is very small one can miss the interesting part regarding the localized impurity state. XES and SXA should be more useful, since they allow for a local projection of the bandstructure by the involvement of a core hole. A La core state can be used to probe the occupied or unoccupied La valence states in XES or SXA, respectively. The signal might be weak because of few La sites, but the spectrum has the advantage of being a selective $`\mathrm{}`$-projection of the local La-DOS. This means that the signal from La should show a narrow band, while a wide unstructured spectrum (similar to the total DOS for the undoped material) is expected if a core level from Sr or an In-impurity is used. Dipolar matrix elements do some reshaping of the pure DOS functions, but they can be estimated quite easily. The question about relaxation around the core hole is more complicated. According to the so-called final-state rule, one observes the spectrum of the relaxed, final state . If so, one should prefer XES, since the relaxed state is the unperturbed ground state, and the complication of a new calculation including a core-hole for SXA, can be avoided.

Theoretical and experimental efforts have been done to determine the degree of band overlap at the X-point in the undoped hexaborides . Calculations using gradient corrections to the LDA or using the GW method, find different values of the band overlap . If the overlap is absent and no doping, the material will be a small-gap semiconductor. The excitonic mechanism can be sensitive to this, and with pressure or measures for changing the $`z`$-parameter one can think of tuning the material for optimal properties. However, the mechanism for magnetism that we have discussed here is not sensitive to the band structure of the undoped case. It is very sensitive to the degree of localization of the La impurity band. If a band calculation gives a too small localization, a wider band with small derivatives of the DOS, then the moment goes to zero. So, unless there is an unknown connection between the band gap for the undoped band structure and the localization of the La impurity band, we believe that the exact shape of the bands in the undoped system is irrelevant for the mechanism for magnetism around a La-impurity. But, since we have seen that the results are sensitive to features of the impurity band, it also is probable that lattice relaxation around the impurity, including variations of the $`z`$-parameter, will have an effect on the magnetic moment.

This work show that the effect of dilute La-doping on the band structure of SrB6 can lead to the appearance of a magnetic moment around the impurity. Unfortunately, temperature variations cannot be studied with the same level of sophistication. The experimental observation of magnetism is interesting, but also its temperature dependence is unususal as magnetic ordering can persist up to about 600 K . The Fermi-Dirac smearing is only one part of the T-depencencies. It is known that thermal structural (vibrational) disorder is important already at 200-300 K. Delicate changes of band structure and the fine details of the DOS due to thermal disorder can be essential for finding a correct description of the physical properties of a material . This, disorder in the distribution of impurities, disorder of the local moments, non-collinearity and orbital magnetism are not considered here, only the electronic excitations given by the Fermi-Dirac distribution. In general, the calculated magnetic moment disappears for $`k_BT`$ in the range 2-4 mRy (300-600 K), but not always. One case with the large cell had the largest moments for quite large T, which could be related to a sharp DOS structures a bit beside the optimal position for EF. If this mean-field picture is sufficient, one would then expect a wide range of ordering temperatures among different samples, like different samples showing different amplitudes of the moments. The effects of spacial irregularities of the impurity distribution of La or other dopants have not been calculated explicitly. But from the separate results with different relative positions between EF and the localized La band, one can imagine that some samples could be apt for large moments at low T, while others with lower moments could resist better at large T.

Finally, the question is whether this mechanism for magnetism is rare and exotic, or if systems other than the hexaborides can show similar effects. From eq. 1 one should have a small total DOS, but large DOS derivative of a sub-band. Pure compounds with wide, hybridized bands can hardly be interesting cases. Dilute doping in a low DOS material seems to give optimal conditions, because $`N_v/N1`$, and an impurity band can localize and become narrow for large distances between impurities. Not all combinations between impurities and host material will do, because the hybridization can make the impurity band to merge into the host bands, as for In in the hexaboride. The requirement of a large $`U_0`$ with the correct sign, will also reduce the number of candidates for magnetism of this kind. Preliminary calculations for metal impurities in some oxides show similar results as for the La doping here, although there is a sizable polarization in the proximity of the impurity .

## IV Conclusion.

We have shown that ferromagnetism occurs in the supercells of LaSr7B48 and LaSr26B162. By doing the calculations with different number of k-points and different basis functions, it has been possible to show that the actual derivative of the local DOS at $`E_F`$, which is necessary for charge transfers, scales with the stability of the magnetic moment. This demonstrates that the mechanism for ferromagnetism is not exclusively found in the exchange energy. The necessary condition of a localized, non-constant La-DOS near $`E_F`$, appears to be enforced as the doping decreases, and the mechanism of charge transfers persists at lower doping levels. Thus, it is likely that the observed weak ferromagnetism for dilute La-doping in hexaborides , is a band effect of the kind that we have described. This conclusion is based on results from first-principle LMTO calculations. The moments are of the order 0.05-0.1 $`\mu _B`$/La for the calculations with the ”best” choice of basis function and the largest number of k-points. However, the fact that the results can be ”boosted” by choosing somewhat extreme linearization energies, brings up the problem of describing a very sensitive band structure with a linearized method. Most band results (total energy, Fermi surface, and the general shape of bands and DOS functions) usually show only a moderate sensitivity to the choice of basis or k-point convergence. But here, although the visual inspection of the DOS show no large differences, there are small differences in the localization of the impurity bands depending on the choice of basis, which turn out to be important for the size of the magnetic moment. Therefore, one must emphasize that the calculations show the likelihood of a ferromagnetic state, but that one can not be affirmative for the absolute values. This situation may reflect the reality, since different samples are reported to have different moments (at most 0.32 $`\mu _B`$/La), while some show only large paramagnetic susceptibilities . Real materials have some randomness in the distribution of the impurities, while the calculations are made for ordered supercells. Clustering and other lattice imperfections can be important for the properties if they lead to a rigid-band like modification of the band structure, or if the impurity states become more localized. From the balance between such effects, it is suggested that optimal conditions for magnetism can be found indirectly by addition of ”rigid-band” dopants like Al. It is probable that other systems, with localized bands from selected dilute doping in low DOS materials, can be found to show similar type of weak magnetism.

###### Acknowledgements.

I am grateful to R. Monnier for discussions, and to H.R. Ott for sending reprints.