text:
# Spin-dependent Bohm trajectories associated with an electronic transition in hydrogen

## 1 Introduction

In the de Broglie-Bohm causal interpretation of quantum mechanics , as opposed to the standard “Copenhagen” interpretation, particles are endowed with well-defined trajectories $`𝐱_i(t)`$ that are determined by the wavefunction $`\psi (𝐱,t)`$ of the quantum system being studied. A compatability with the statistical results of quantum mechanics is achieved by assigning an uncertainty in the initial conditions of the particles according to the probability density function $`\rho (x,0)=|\psi (x,0)|^2`$.

The causal interpretation continues to receive attention as a way of addressing the question of the incompleteness of standard quantum mechanics (see, for example, ). Recently, however, causal trajectories have been playing a more significant role in practical calculations of chemical physics and quantum chemistry, for example (to name only a few references), quantum tunnelling dynamics , nonadiabatic transitions , reactive scattering , dissociation dynamics  and hybrid classical/quantum schemes to study complex systems .

Bohm himself  introduced the idea of studying transitions in terms of the causal interpretation, examining the Franck-Hertz experiment and the photoelectric and Compton effects. He attempted to show that the seemingly discontinuous and poorly defined transfers of energy and momentum in transitions could be accounted for in a continuous matter by means of the “quantum potential” that arises in the causal formalism. More recently in this vein, Dewdney and Lam  studied transitions of (spinless) particles in a one-dimensional infinite square well potential. In this paper we employ the de Broglie-Bohm deterministic approach to study the problem of a $`1s`$-$`2p`$ electronic transition in hydrogen induced by an oscillating (semiclassical) electromagnetic field, taking the spin of the electron into account.

It is instructive to review briefly the main ideas of the causal interpretation. First, the wavefunction $`\psi `$ for a particle is written in the form

$$\psi (𝐱,t)=R(𝐱,t)e^{iS(𝐱,t)/\mathrm{}},$$
(1)
where $`R`$ and $`S`$ are real-valued. Substitution of Eq. (1) into the time-dependent Schrödinger equation,

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}^2\psi +V\psi ,$$
(2)
yields the following coupled equations in $`R`$ and $`S`$:

$$\frac{R^2}{t}+\left(R^2\frac{S}{m}\right)=0$$
(3)
and

$$\frac{S}{t}+\frac{(S)^2}{2m}+V+Q=0,$$
(4)
where

$$Q=\frac{\mathrm{}^2}{2m}\frac{^2R}{R}$$
(5)
is called the quantum potential. Eq. (3) is the standard continuity equation of quantum mechanics. It can be viewed as governing the evolution of a compressible, irrotational fluid with density $`\rho =\psi ^{}\psi =R^2`$ and velocity $`𝐯=S/m`$ as, indeed, was done by Madelung . Eq. (4) has the same form as the classical Hamilton-Jacobi equation for a particle that moves under the influence of potentials $`V`$ and $`Q`$. Bohm’s unique interpretation of these equations was that the particle has a well-defined trajectory defined by the quantum equation of motion

$$𝐩=m\dot{𝐱}=S.$$
(6)
As stated earlier, compatability with standard quantum mechanics is achieved by viewing the initial conditions of trajectories as “hidden variables” with associated uncertainties as described by the probability distribution $`\rho (𝐱,0)`$.

Recall that the Schrödinger current j associated with the wavefunction $`\psi `$ is given by

$$𝐣=\frac{\mathrm{}}{2mi}[\psi ^{}\psi \psi \psi ^{}].$$
(7)
A comparison of Eqs. (6) and (7) shows that

$$𝐣=\frac{1}{m}\rho 𝐩.$$
(8)
The momentum defined in Eq. (6) is not unique in generating the same statistical predictions as quantum mechanics. Holland  has shown that Eqs. (6) and (8) apply only to spinless particles. For particles with spin, an additional term is necessary in order that the Schrödinger equation of motion be consistent with a relativistic formulation. The condition of Lorentz invariance implies that the momentum of a particle with spin s, even in the non-relativistic limit, must be given by

$$𝐩=S+\mathrm{log}\rho \times 𝐬.$$
(9)
The associated current

$$𝐣=\frac{1}{m}\rho S+\frac{1}{m}\rho \times 𝐬,$$
(10)
has been referred to as the Pauli current, the nonrelativistic limit of the Dirac current, as opposed to Eq. (8) which is the nonrelativistic limit of the Gordon current . Consistency with Dirac theory requires that the Schrödinger equation be regarded as describing an electron in an eigenstate of spin . As regards the causal interpretation, the spin-dependent term was also discussed in  but only for the Pauli equation and not the Schrödinger equation.

In the case of a hydrogen atom, the momentum equation (6) predicts that an electron in any real eigenstate will be stationary since $`S=0`$. This counterintuitive result no longer applies when Eq. (9) is used. For example, consider an electron in the $`1s`$ ground eigenstate with wavefunction

$$\psi _{100}=\frac{1}{\sqrt{\pi a^3}}e^{r/a},$$
(11)
where $`a=\mathrm{}^2/(me^2)`$ is the Bohr radius. Also assume that the electron is in a definite spin eigenstate so that its spin vector is given by $`𝐬=\frac{\mathrm{}}{2}𝐤`$. Then the existence of the term $`\mathrm{log}\rho \times 𝐬`$ in Eq. (9) implies that the electron’s polar coordinates $`r`$ and $`\theta `$ are constant and that the angle $`\varphi `$ evolves in time as

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{mar}.$$
(12)
Therefore, the electron revolves about the $`z`$-axis at a constant angular velocity.

We also state, for future reference, the result for an electron in the (real) $`2p_0`$ eigenstate

$$\psi _{210}=\frac{1}{\sqrt{32a^5}}re^{r/2a}\mathrm{cos}\theta ,$$
(13)
again with spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$: The polar coordinates $`r`$ and $`\theta `$ are again constant and the angle $`\varphi `$ evolves as

$$\frac{d\varphi }{dt}=\frac{\mathrm{}}{2mar}.$$
(14)
Note that the angular velocity is one-half that of a $`1s`$ ground state electron. In  we examined the trajectories of electrons in a number of hydrogenic eigenstates. From Eq. (9), these trajectories must lie on level curves of both $`|\psi |`$ and $`z`$ and revolve about the $`z`$-axis with constant angular velocity.

In this paper we examine solutions of the equation of motion (9) for an electron with spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$ (the “$`\alpha `$” or “spin up” state) as it undergoes a transition from the $`1s`$ to $`2p_0`$ state in hydrogen due to the presence of an oscillating electric field. We may assume that the electron has constant spin vector since the hamiltonian describing the atom in the field (see next section) is spin-independent. The wavefunction of the electron $`\mathrm{\Psi }(𝐱,𝐬,t)`$ may then be written as the tensor product $`\psi (𝐱,t)\zeta (𝐬)`$ where $`\zeta (𝐬)`$ is assumed to be an eigenfunction of the commuting spin operators $`\widehat{S}^2`$ and $`\widehat{S}_z`$, with $`\widehat{S}^2\zeta =\frac{3\mathrm{}}{4}\zeta `$ and $`\widehat{S}_z\zeta =\frac{\mathrm{}}{2}\zeta `$. As such, the remainder of our discussion may simply be focussed on the evolution of the spatial portion of the wavefunction, $`\psi (𝐱,t)`$ according to Eq. (2).

As in the case of the examples listed above, the momentum term $`\mathrm{log}\rho \times 𝐬`$ will be responsible for the revolution of the electron about the $`z`$-axis. This is accompanied by a complicated motion over a hyperboloid surface of revolution that is determined from the functional forms of the $`1s`$ and $`2p_0`$ wavefunctions. Moreover, the course of the transition from the ground state to the excited state can be characterized by looking at the energy of the electron along the trajectory and the angular velocity of the revolution about the $`z`$-axis. The energy and $`\varphi `$-angular velocity of the electron evolve from $`1s`$ ground state values to $`2p`$ excited state values.

In , as a precursor to this study, we examined the trajectories dictated by Eq. (9) for an electron with spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$ and spatial wavefunction that begins as a linear combination of $`1s`$ and $`2p`$ hydrogenic eigenfunctions:

$$\psi (𝐱,0)=c_1\psi _{100}(𝐱)+c_2\psi _{210}(𝐱),$$
(15)
where $`|c_1|^2+|c_2|^2=1`$. The time evolution of this wavefunction under the hydrogen atom hamiltonian is

$$\psi (𝐱,t)=c_1\psi _{100}(𝐱)e^{iE_1t/\mathrm{}}+c_2\psi _{210}(𝐱)e^{iE_2t/\mathrm{}}.$$
(16)
Many of the qualitative features of the $`1s`$-$`2p_0`$ transition problem studied below are captured by this model, most notably the invariant hyperboloid surfaces of revolution on which trajectories lie. As expected, however, the more detailed time evolution of the electron trajectories over these surfaces due to the oscillating field is missing.

## 2 Solution of the transition problem

The hamiltonian used to describe this transition will have the form $`\widehat{H}=\widehat{H}_0+\widehat{H}^{}`$, where $`\widehat{H}_0`$ is the hydrogen atom hamiltonian and $`\widehat{H}^{}`$ represents an oscillating electric field $`𝐄=E_0\mathrm{cos}\omega t𝐤`$. It can be shown  that if $`\omega `$ is chosen to be sufficiently close to the $`1s`$-$`2p`$ transition frequency

$$\omega _0=\frac{E_2E_1}{\mathrm{}},$$
(17)
so that $`\omega _0\omega \omega _0+\omega `$, then the hamiltonian representing the semiclassical radiation, $`\widehat{H}^{}=qzE_0\mathrm{cos}\omega t`$, is well approximated by

$$\widehat{H}^{}=\frac{1}{2}qzE_0e^{i\omega t},$$
(18)
since the term $`i(\omega +\omega _0)^1\mathrm{sin}\omega t`$ is negligible. Here, $`q`$ denotes the electric charge and $`\omega _01.549\times 10^{16}\text{s}^1`$. This approach allows the equations for the wavefunction coefficients to be solved exactly so that perturbation methods need not be employed. The closeness of $`\omega `$ to $`\omega _0`$ also allows the transition probability to approach unity at various times rather than remaining small for all times.

The probability of transition between two states $`|\psi _1`$ and $`|\psi _2`$ is related to the matrix element $`\psi _2\left|\widehat{H}^{}\right|\psi _1`$. In the case of the hydrogen atom the only nonvanishing matrix element is between the ground state $`\psi _{100}`$ and the $`2p_0`$ state $`\psi _{210}`$:

$`\psi _{100}\left|\widehat{H}^{}\right|\psi _{210}`$ $`=`$ $`\psi _{100}\left|qE_0r\mathrm{cos}\theta \right|\psi _{210}{\displaystyle \frac{1}{2}}e^{i\omega t}`$ (19)
$`=`$ $`{\displaystyle \frac{64\sqrt{2}}{243}}aqE_0e^{i\omega t}.`$
The time-dependent wavefunction can be written as a linear combination of $`|\psi _{100}`$ and $`|\psi _{210}`$:

$$\psi (t)=c_a(t)\psi _{100}e^{iE_1t/\mathrm{}}+c_b(t)\psi _{210}e^{iE_2t/\mathrm{}}.$$
(20)
Substitution into the Schrödinger equation yields the following equations for $`c_a(t)`$ and $`c_b(t)`$:

$$\begin{array}{cc}\hfill \dot{c_a}& =\frac{i}{\mathrm{}}\frac{V_{12}}{2}e^{i(\omega _0\omega )t}c_b,\hfill \\ \hfill \dot{c_b}& =\frac{i}{\mathrm{}}\frac{V_{12}}{2}e^{i(\omega _0\omega )t}c_a,\hfill \end{array}$$
(21)
where

$$V_{12}=\frac{128\sqrt{2}}{243}aqE_0.$$

These DEs can be solved exactly to give

$$\begin{array}{cc}\hfill c_a(t)& =\frac{\sigma +\mathrm{\Omega }}{2\sigma }e^{\frac{1}{2}i(\mathrm{\Omega }\sigma )t}+\frac{\sigma \mathrm{\Omega }}{2\sigma }e^{\frac{1}{2}i(\mathrm{\Omega }+\sigma )t},\hfill \\ \hfill c_b(t)& =\frac{\nu }{2\sigma }e^{\frac{1}{2}i(\mathrm{\Omega }\sigma )t}\frac{\nu }{2\sigma }e^{\frac{1}{2}i(\mathrm{\Omega }+\sigma )t},\hfill \end{array}$$
(22)
where

$$\begin{array}{cc}\hfill \mathrm{\Omega }& =\omega _0\omega ,\hfill \\ \hfill \nu & =\frac{V_{12}}{\mathrm{}},\hfill \\ \hfill \sigma & =\sqrt{\mathrm{\Omega }^2+\nu ^2}.\hfill \end{array}$$
(23)
The wavefunction $`\psi (𝐱,t)`$ may now be written explicitly as

$$\psi (𝐱,t)=\frac{1}{\sqrt{\pi a^3}}c_a(t)e^{r/a}e^{iE_1t/\mathrm{}}+\frac{1}{\sqrt{32\pi a^5}}c_b(t)re^{r/2a}\mathrm{cos}\theta e^{iE_2t/\mathrm{}}.$$
(24)
To compute the momentum according to Eq. (9), note that the wave function, and hence first term $`S`$, has only $`\widehat{r}`$ and $`\widehat{\theta }`$ components. Since we are assuming a constant spin vector $`𝐬=\frac{\mathrm{}}{2}\widehat{k}`$, it follows that the vector $`\mathrm{log}\rho \times 𝐬`$ points in the $`\widehat{\varphi }`$ direction.

Calculating $`S`$ from Eq. (24) yields

$$\begin{array}{cc}\hfill p_{\widehat{r}}& =\frac{\mathrm{}\nu \beta }{2\sigma }\frac{\mathrm{cos}\theta e^{3r/2a}(1+\frac{r}{2a})T(t)}{D(r,\theta ,t)},\hfill \\ \hfill p_{\widehat{\theta }}& =\frac{\mathrm{}\nu \beta }{2\sigma }\frac{\mathrm{sin}\theta e^{3r/2a}T(t)}{D(r,\theta ,t)},\hfill \end{array}$$
(25)
where $`\beta =4\sqrt{2}a`$ is the ratio of the normalizing factors of the two wavefunctions,

$$T(t)=\mathrm{cos}\omega _0t\mathrm{sin}\sigma t\frac{\mathrm{\Omega }}{\sigma }\mathrm{sin}\omega _0t+\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\sigma t\mathrm{sin}\omega _0t$$
(26)
and

$$\begin{array}{c}D(r,\theta ,t)=e^{2r/a}\frac{1}{2\sigma ^2}(\sigma ^2+\mathrm{\Omega }^2+\nu ^2\mathrm{cos}\sigma t)+\beta ^2r^2e^{r/a}\mathrm{cos}^2\theta \frac{\nu ^2}{2\sigma ^2}(1\mathrm{cos}\sigma t)\hfill \\ \hfill +\frac{\nu }{\sigma }\beta re^{\frac{3r}{2a}}\mathrm{cos}\theta (\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\omega _0t\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\sigma t\mathrm{cos}\omega _0t\mathrm{sin}\omega _0t\mathrm{sin}\sigma t).\end{array}$$
(27)
The denominator $`D(r,\theta ,t)`$ in the above expressions is proportional to $`|\psi (𝐱,t)|^2`$.

The second term in Eq. (9), $`\mathrm{log}\rho \times 𝐬`$, can be computed in the (right-handed) spherical polar coordinate system

$$𝐀\times 𝐁=\left|\begin{array}{ccc}\widehat{\theta }& \widehat{\varphi }& \widehat{r}\\ A_\theta & A_\varphi & A_r\\ B_\theta & B_\varphi & B_r\end{array}\right|.$$
(28)
Using the simplification

$$\mathrm{log}\rho =2\text{Re}\left(\frac{(\psi )\psi ^{}}{\psi ^{}\psi }\right),$$

we find the $`\widehat{\varphi }`$ component of the momentum to be

$$p_{\widehat{\varphi }}=\frac{\mathrm{}\beta }{D}(\chi _r\mathrm{sin}\theta \chi _\theta \mathrm{cos}\theta ),$$
(29)
where

$`\chi _r`$ $`={\displaystyle \frac{1}{\beta a}}|c_a|^2e^{2r/a}+\beta |c_b|^2\mathrm{cos}^2\theta e^{r/a}r(1{\displaystyle \frac{r}{2a}})+\mathrm{cos}\theta e^{3r/2a}(1{\displaystyle \frac{3r}{2a}})T^{}(t),`$
$`\chi _\theta `$ $`=\beta |c_b|^2e^{r/a}\mathrm{sin}\theta \mathrm{cos}\theta re^{3r/2a}\mathrm{sin}\theta T^{}(t)`$
and

$$T^{}(t)=\frac{\nu }{2\sigma }(\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\omega _0t\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\sigma t\mathrm{cos}\omega _0t\mathrm{sin}\sigma t\mathrm{sin}\omega _ot).$$
(30)
In summary, the three components of the momentum are given by Eqs. (25) and (29). It is worth noting that the spin-dependent momentum term $`\mathrm{log}\rho \times 𝐬`$ in Eq. (9) is responsible for the $`\varphi `$-momentum $`p_{\widehat{\varphi }}`$.

It is useful to rescale these equations by defining the following dimensionless variables:

$$\xi =\frac{r}{a},\tau =\omega _0t.$$
(31)
In these variables, Eqs. (25) and (29) give rise to the following system of differential equations in $`\xi `$, $`\theta `$, and $`\varphi `$ as functions of $`\tau `$:

$$\begin{array}{cc}\hfill \frac{d\xi }{d\tau }& =\frac{\nu }{3\sqrt{2}\sigma }\left(\mathrm{cos}\theta e^{3\xi /2}(1+\frac{\xi }{2})\right)\frac{\stackrel{~}{T}(\tau )}{D}\hfill \\ \hfill \frac{d\theta }{d\tau }& =\frac{\nu }{3\sqrt{2}\sigma }\left(\mathrm{sin}\theta \frac{e^{3\xi /2}}{\xi }\right)\frac{\stackrel{~}{T}(\tau )}{D}\hfill \\ \hfill \frac{d\varphi }{d\tau }& =\frac{\nu }{3\sqrt{2}\sigma \xi D}(\chi _r+\chi _\theta \mathrm{cot}\theta ),\hfill \end{array}$$
(32)
where we have used the following relations:

$$\frac{d\varphi }{d\tau }=\frac{p_\varphi }{ma\omega _0\xi \mathrm{sin}\theta }$$

and

$$\frac{\mathrm{}}{ma\omega _0}=\frac{8}{3}a.$$

Some important qualitative features of the solutions to these DEs may be extracted. First note that the $`\xi `$ and $`\theta `$ DEs are independent of $`\varphi `$. ¿From these two equations, we have

$$\frac{d\xi }{d\theta }=\xi \left(1+\frac{\xi }{2}\right)\mathrm{cot}\theta .$$
(33)
This DE is easily solved to give

$$\xi =\frac{2}{A\mathrm{sin}\theta 1},A=\frac{2+\xi _0}{\xi _0\mathrm{sin}\theta _0}>1,$$
(34)
where $`\xi _0=\xi (0)`$ and $`\theta _0=\theta (0)`$. This relation defines a family of hyperbolae in vertical planes that contain the $`z`$-axis. i.e., $`\varphi =\varphi _0`$, a constant. (Note that the right hand side of Eq. (33) is determined by the functional forms of the $`1s`$ and $`2p_0`$ wavefunctions. Other allowable pairs of wavefunctions will yield different types of curves.) The solutions of Eq. (32) therefore lie on surfaces that are obtained by rotating these hyperbolae about the $`z`$-axis (see  for a more detailed discussion). However, the time-dependent behaviour of the trajectories lying on these invariant sets must be determined numerically.

## 3 Numerical results

In choosing the parameters for the numerical integration, there are various factors that must be taken into account. First, recall that in order to use the hamiltonian $`\widehat{H}^{}=\frac{1}{2}qzE_0e^{i\omega t}`$, we require that the perturbing frequency $`\omega `$ be close to the transition frequency $`\omega _0`$ so that $`\mathrm{\Omega }=\omega _0\omega \omega _0+\omega `$. Therefore, we cannot allow $`\mathrm{\Omega }`$ to be too large, i.e., $`\mathrm{\Omega }O(10^{13})`$.

If we wish to be fairly confident that a transition will occur, it is also necessary that the coefficient $`c_b`$ become large in magnitude at some time, i.e., roughly unity. Recall from Eq. (22) that $`c_b\nu /2\sigma `$, and $`\nu E_0`$. Also, the derivatives in Eq. (32) are proportional to $`\nu /2\sigma `$, and it is desirable for numerical integration that they evaluate to order 1. Therefore, $`|\nu |`$ should be significant compared to $`2\sigma `$. Now, recall that $`\sigma =\sqrt{\mathrm{\Omega }^2+\nu ^2}`$ so that

$$\frac{|\nu |}{2\sigma }=\frac{1}{2}\left[\frac{\mathrm{\Omega }^2}{\nu ^2}+1\right]^{1/2}.$$

Therefore we require that $`|\nu |`$ not be too large. Because $`|\nu |`$ is proportional to the field strength, we are free to choose a small value.

Another factor to consider is that there are two angular frequencies in the problem, namely $`\sigma `$ and $`\omega _0`$ from Eqs. (26) and (30). For numerical stability it is best if these are within several orders of magnitude of each other. Therefore, $`\sigma `$ should not be too small in comparison with $`\omega _0`$. It is also desirable that $`c_b`$ become significant in a reasonable time. From Eq. (22), we find that

$$|c_b(t)|^2=\left(\frac{\nu }{\sigma }\right)^2\mathrm{sin}^2\frac{\sigma }{2}t.$$
(35)
Therefore, having $`\sigma `$ appropriately scaled will mean that there is a high probability of seeing a transition within a reasonable time.

The above considerations imply that:

1. $`|\nu |`$ cannot be too large because $`[\mathrm{\Omega }^2/\nu ^2+1]^{1/2}`$ must be O(1),
2. $`\sigma =\sqrt{\mathrm{\Omega }^2+\nu ^2}`$ cannot be more than several orders of magnitude smaller than $`\omega _o`$ and
3. we cannot increase $`\sigma `$ by increasing $`\mathrm{\Omega }`$, because we require that $`\mathrm{\Omega }O(10^{13})`$.

With these points in mind, the parameters have been chosen as follows:

$$\begin{array}{cc}\hfill E_0& =8.8\times 10^7\text{V/m},\hfill \\ \hfill \mathrm{\Omega }& =1.55\times 10^{12}\text{s}^1,\hfill \end{array}$$
(36)
so that

$$\begin{array}{cc}\hfill \nu & =5.1\times 10^{12},\hfill \\ \hfill \sigma & =5.32\times 10^{12}.\hfill \end{array}$$
(37)
The results of numerical integration for two choices of initial conditions are shown in Figs. (1) and (2). Figs. (3) and (4) show the same trajectories split over five consecutive time intervals. In these latter plots, the hyperboloid surfaces of revolution on which the trajectories lie, described at the end of the previous section, are discernable.

In standard semiclassical treatments of this problem (see, for example, , pp. 282-285), the maximum probability of transition occurs when $`\mathrm{sin}^2\frac{1}{2}\mathrm{\Omega }t=1`$, or $`\mathrm{\Omega }t=(2k+1)\pi `$ for $`k=0,1,2,\mathrm{}`$. In our formulation, from Eq. (35) the maximum probability of transition occurs when $`\mathrm{sin}^2\frac{1}{2}\sigma t=1`$, or $`\sigma t=(2k+1)\pi `$ for $`k=0,1,2,\mathrm{}`$. In other words, the dependence on $`\mathrm{\Omega }`$ is replaced with one on $`\sigma =\sqrt{\mathrm{\Omega }^2+\nu ^2}`$. This is due to our choice of hamiltonian and the fact that the equations in (22) result from an exact integration of Eq. (21) rather than a perturbation approach. ¿From the above, the first occurrence of the maximum probability of transition in our problem will occur at $`t=\pi /\sigma `$ or, in dimensionless time, $`\tau =\pi \omega _0/\sigma `$. ¿From the numerical values chosen for the parameters above, this corresponds to $`\tau 9000`$.

The causal interpretation allows us to examine the process in more detail by looking at the angular velocity $`\dot{\varphi }`$ exhibited by the trajectories. At time $`t=\tau =0`$, each trajectory begins at the $`1s`$ ground state wavefunction, viz., $`c_a(0)=1`$ and $`c_b(0)=0`$, and we expect the electron to revolve about the $`z`$-axis with angular frequency given by Eq. (12). When the wavefunction is approximately equal to the $`2p_0`$ eigenstate, viz., $`c_a0`$ and $`c_b1`$, then we expect $`\dot{\varphi }`$ to be given roughly by Eq. (14).

Since the computations were performed in scaled variables $`\xi `$ and $`\tau `$, cf. Eq. (31), we must first rescale the $`1s`$ and $`2p_0`$ angular velocities in $`\varphi `$, cf. Eqs. (12) and (14), in order to assess the numerical results. For the $`1s`$ state, Eq. (12) becomes

$$\frac{d\varphi }{d\tau }=\frac{\mathrm{}}{ma^2\omega _0\xi }=\frac{8}{3\xi }.$$
(38)
The angular velocity for the $`2p_0`$ state is simply one-half the above result.

The values of the scaled angular velocity $`d\varphi /d\tau `$ associated with the trajectory shown in Fig. (1) are presented as a function of $`\tau `$ in Fig. (5). (In Fig. (6) are shown the corresponding values of $`\varphi `$.) Recall that this trajectory began with the value $`\xi (0)=4`$. ¿From Fig. (1), $`d\varphi /d\tau `$ is observed to be roughly $`2/3`$ near $`\tau =0`$. This is in agreement with Eq. (38) – the electron is revolving about the $`z`$-axis at the $`1s`$ rate.

Near $`\tau =9000`$, roughly the time for $`|c_b(\tau )|^2`$ to achieve its maximum value, we observe that $`d\varphi /d\tau 0.3`$. At that time, $`\xi 4.5`$. This is in agreement with the scaled $`2p_0`$ rate $`4/(3\xi )0.3`$.

In the causal interpretation it is also possible to ascribe a value of energy along a trajectory and thereby examine the time-dependent behaviour of the energy during a transition. Note that this is in contrast to standard quantum mechanics, in which we can only compute probabilities of measuring energy eigenstates $`E_1`$ or $`E_2`$ after the perturbation has been turned off. One method of computing the energy is to use the result $`E=S/t`$ implied by the Hamilton-Jacobi equation (4). However, even when applied to the relatively simple wavefunction $`\psi `$ in Eq. (24), this method is quite cumbersome. On the other hand, Holland’s method of assigning values of an observable quantity $`A`$ to points on a trajectory, which we outline below (see , pp. 91-93), is quite easily implemented in our problem.

The usual expectation value of an observable $`\widehat{A}`$ is

$$<\widehat{A}>=\psi \left|\widehat{A}\right|\psi =\frac{\psi ^{}\widehat{A}\psi d^3𝐱}{\psi ^{}\psi d^3𝐱}.$$

Since $`\widehat{A}`$ must be Hermitian, ultimately only the real portion of the expression contributes to the integral. Hence define the local expectation value to be

$$A(𝐱,t)=\frac{\text{Re}\{\psi ^{}\widehat{A}\psi \}}{\psi ^{}\psi }.$$
(39)
To find the local energy expectation value along the trajectory, we compute this quantity using the hamiltonian $`\widehat{H}=\widehat{H}_0+\widehat{H}^{}`$ where, as before, $`\widehat{H}_0`$ is the usual hydrogen hamiltonian, and

$$\widehat{H}^{}=\frac{1}{2}eE_0e^{i\omega t}(r\mathrm{cos}\theta ).$$

Thus

$$E(𝐱,t)=\frac{\text{Re}\{\psi ^{}(\widehat{H}_0+\widehat{H}^{})\psi \}}{\psi ^{}\psi }=\frac{\text{Re}\{\psi ^{}\widehat{H}_0\psi \}}{\psi ^{}\psi }+\frac{\text{Re}\{\psi ^{}\widehat{H}^{}\psi \}}{\psi ^{}\psi },$$
(40)
where $`\psi (𝐱,t)`$ is given by Eq. (20). Note that the relationship to the classical energy is not direct; this local energy expectation value is not a first integral of motion, simply a quantity whose average, over trajectories, yields the quantum mechanical energy expectation value. Here we will refer to it as the local energy. After some manipulations, the first term in Eq. (40) is given by

$$\frac{\text{Re}\{\psi ^{}(\widehat{H}_0)\psi \}}{\psi ^{}\psi }=\frac{|c_a|^2E_1\psi _{100}^2+|c_b|^2E_2\psi _{210}^2+\psi _{100}\psi _{210}(E_2+E_1)\text{Re}\{c_ac_be^{i\omega t}\}}{\psi ^{}\psi }$$

and the second term is simply

$$\frac{\text{Re}\{\psi ^{}\widehat{H}^{}\psi \}}{\psi ^{}\psi }=\frac{\text{Re}\{\psi ^{}(\frac{1}{2}eE_0e^{i\omega t}(r\mathrm{cos}\theta )\psi )\}}{\psi ^{}\psi }=\frac{1}{2}eE_0r\mathrm{cos}\theta \mathrm{cos}\omega t$$

where, referring to Eq. (30),

$$\text{Re}\{c_ac_be^{i\omega t}\}=\frac{\nu }{2\sigma }(\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\omega _0t\frac{\mathrm{\Omega }}{\sigma }\mathrm{cos}\sigma t\mathrm{cos}\omega _0t\mathrm{sin}\sigma t\mathrm{sin}\omega _ot)=T^{}(t).$$

Therefore we have

$$E=\frac{|c_a|^2E_1\psi _{100}^2+|c_b|^2E_2\psi _{210}^2+\psi _{100}\psi _{210}(E_2+E_1)T^{}(t)}{\psi ^{}\psi }\frac{1}{2}eE_0r\mathrm{cos}\theta \mathrm{cos}\omega t.$$
(41)
This function can be evaluated along the trajectories shown in Figs. (1) and (2). The results for the first trajectory are shown in Fig. (7).

Note that at $`\tau =0`$, the local energy is equal to the $`1s`$ ground state energy $`13.6\text{eV}`$, and that after the transition, the energy oscillates near the value of $`E_2=3.4\text{eV}`$, corresponding to the energy of hydrogenic $`n=2`$ states. The energy is not constant even in a proximity of these points because of the extra oscillating perturbation $`\widehat{H}^{}`$ that represents the semi-classical radiation field. The spike in the energy corresponds to a point along the trajectory that approaches a zero of the wavefunction. Intuitively, one can understand the appearance of such a spike because the electron cannot spend time in regions where $`\psi ^{}\psi `$ is very small – it receives a ‘kick’ from the quantum potential which corresponds to a sudden rise in its energy. Plots of the quantum potential for several systems and a discussion are given in Holland  and references within, for example, . The energy along the trajectory shown in Fig. (2) is shown in Fig. (8).

Finally, we mention that we have computed trajectories to higher times. At $`\tau 18000`$, the scaled angular velocity $`d\varphi /d\tau `$ is observed to oscillate about the $`1s`$ value. The energy is also observed to oscillate about the $`1s`$ value of -13.6 eV. This is in accordance with $`|c_b(\tau )|^20`$ from Eq. (35) and the values of the parameters used in the calculation. We conclude that the transition has reversed and that the electron has returned to the $`1s`$ ground state.

## 4 Concluding remarks

We have studied the problem of a simple $`1s`$ to $`2p_0`$ electronic transition in hydrogen – induced by an oscillating (semiclassical) radiation field – in terms of the causal interpretation of quantum mechanics. In contrast to Bohm’s original formulation $`𝐩=S`$, however, we have employed an additional spin-dependent term $`\mathrm{log}\rho \times 𝐬`$ in the equation of motion, where $`\rho =\psi ^{}\psi `$. The electron is assumed to be in a “spin up” eigenstate, with associated spin vector $`𝐬=\frac{\mathrm{}}{2}𝐤`$, during the transition. The electronic trajectories lie on invariant hyperboloid surfaces of revolution about the $`z`$-axis. The nature of the trajectories over these surfaces is quite complex due to the quasiperiodic nature of the equations of motion which, in turn, arise from the interplay of the $`1s`$-$`2p`$ transition frequency $`\omega _0`$, the frequency $`\omega `$ of the oscillating radiation and $`E_0`$, the field strength.

As the electron moves over the invariant surface it also revolves about the $`z`$-axis due to the spin-dependent momentum term. The progress of the transition can be tracked by observing the local energy $`E`$ and the angular velocity $`\dot{\varphi }`$ of the electron. Beginning at values associated with the $`1s`$ ground state, both quantities are seen to evolve toward values associated with the $`2p_0`$ excited state in time. The causal interpretation has offered a way to examine the phenomenon of transition which is not limited to computing the probabilities of measuring $`E_1`$ and $`E_2`$ at various times after the perturbation has been turned off.

## Acknowledgments

We gratefully acknowledge that this research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) in the form of a Postgraduate Scholarship (CC) and an Individual Research Grant (ERV). CC also acknowledges partial financial support from the Province of Ontario (Graduate Scholarship) as well as the Faculty of Mathematics, University of Waterloo.


format: text

text:
# Generalised Modal Analysis with the Padé-Laplace transform

## I Introduction

Decomposing a series of data points depending on a single parameter (that might be time, temperature, pressure…) is an important problem in Science and Technology.

The data might be conductivity, chemical concentration, stock market volume, a communication signal etc… It is required that the data is suspected to contain an unknown number of exponential components possessing complex amplitudes and having complex arguments linear functions of the single parameter at hand tufts ; cadzow ; rao .

The important point is that we do not have to specify the number of exponential components: the Padé-Laplace algorithm finds it. This means, the algorithm performs Fourier analysis in the case of periodic or multiperiodic input data, decomposes a purely decaying function into a set of exponentials having real arguments in contrast to Fourier analysis where the arguments are pure imaginary, and finally analyses a function that decays while simultaneously oscillating as consisting of a set of exponentials having complex arguments. In all three cases the amplitudes are determined in complex form. It is a well-known fact that Fourier/spectral analysis breaks down when dealing with exponentially decaying functions and numerical analysis or standard fitting techniques fail when dealing with functions of that sort, because of the ill-posed nature of the problem of exponential fitting acton . G-MAPLE circumvents ill-posedness and succesfully gives, not only the number of exponential components in the function, but also the corresponding amplitudes and arguments.

The applications of Generalised Modal analysis are numerous. Some of the fields of applications are:

1. Numerical fitting software:
Exponential and trigonometric fitting of functions as used in the general Pade rational approximation approach encountered in Quantum Mechanics leung , Particle physics, statistical physics wittle ….
2. Time series software:
Determination of time constants in arbitrary profiles of use in general diagnosis, seasonal trends, forecasting etc…scharf .
3. Communications engineering:
Separation of a signal from its echoes in telephone lines, multipath links (line of sight microwave propagation or mobile outdoor/ indoor); traffic analysis requiring the determination of different parameters yielding traffic volume such as mean arrival times at a given channel sarkar .
4. Chemical engineering, oil, gas, pharmaceutical…:
Determination of the order of a given chemical reaction, detection of intermediate reactions, help in the finding of intermediate reactants or processes etc… here, one expects exponential functions of temperature as well, due to the various thermal activation or propagation processes rasol .
5. Control systems:
Detection of the optimal or hidden number of parameters governing the behaviour of a complex system or plant…
6. Mechanical engineering:
Determination of the number of degrees of freedom (linear or angular displacements or velocities) in the the response of mechanical system such as a robot, mechanical arm, calculation of friction constants controlling the motion etc…
7. Population models in economics, biology, ecology…:
Determination of number and type of species (biological or economic entities) and understanding the underlying dynamics giving rise to the overall observed population or stock market volume at any given time… For instance, it is used in fluorescence decay studies in Biochemistry foucault ; scalley .
8. Speech and acoustic engineering:
Determination of basic components of a phoneme in the time domain, in contrast to the windowed spectral analysis method that is usually performed. echo and reverberation investigation in concert halls…
9. Geophysical prospecting and seismic data processing:
Separation of the various echoes produced by a shock wave, seism and identification of the various types and spatial extent of different layers of materials, strata… mehta .
10. Micro-electronics industry:
Determination of types of impurities or dopants when one monitors the total current versus temperature or the time dependent capacitance in DLTS (Deep Level Transient Spectroscopy) studies…

More recently, it was applied to gene sequencing where exponential functions are used in the calculation of partition functions of DNA loop sequences gorban ; yeram00 .

## II Software presentation

G-MAPLE is a software that decomposes the given data in a set of exponential functions having complex amplitudes and arguments. The novelty is that G-MAPLE determines the unknown number of exponentials in the data along with their complex amplitudes and arguments yeram87 .
G-MAPLE does a Taylor expansion of the Laplace transform of the data around a selected point in the complex plane, fits the expansion to a rational approximation, finds its poles and infers from them the arguments of the exponential components and the values of their corresponding amplitudes. The number of exponential components is found by increasing steadily the degree of the rational approximation. Presumably a set of poles will appear beyond a certain degree and keep appearing later on, regardless of the degree. The number of stable poles is the number of exponential components sought. These stable poles and their associated amplitudes give the exponential functions the sum of which is the best fit to the initial data.

Hence the plan of the G-MAPLE strategy claverie  is:

* Preprocessing of the data: Filtering/Estimating the noise contained in the data or/and altering the sampling rate of the data.
* Determination of the optimal point $`p_0`$ for the Taylor expansion of the Laplace transform of the data.
* Padé rational approximation to the Taylor expansion.
* Inspection of pole stability.
* Consistency analysis.

### II.1 Preprocessing of the data:

This is performed with a spline program. Once the data is entered an an estimate of the noise and a smoothing parameter are provided. Altering and trying several sampling rates is done because it is important for the accuracy required in the calculation of the Taylor coefficients.

### II.2 Determination of the optimal expansion point $`p_0`$:

1. Initial study of the approximate locations of the poles from the Laplace transform of the data.
2. Performing the Taylor expansion for a given degree around several expansion points and examining the behaviour of the coefficients with the degree.

The estimation of the Taylor coefficients of the Laplace transform claverie  of the data is done for several trial p value. An acceptable value for the expansion point $`p_0`$ should yield a set of coefficients whose magnitude decrease smoothly with their algebraic normalised value staying at about the same level. This is done by a simple inspection of the plots of the above sets of coefficients versus the various expansion points. Once the point $`p_0`$ is found one can start the rational approximation analysis. By the same token the damping coefficient is calculated as well. This is required for rapidly varying data. It is possible sometimes to alter slightly the damping coefficient in order to get a better expansion point.

### II.3 Rational approximation procedure:

Once the Taylor coefficients of the Laplace transform have been determined about the optimal expansion point $`p_0`$, one has to vary the degree of the rational approximation to the Taylor expansion. One selects a high degree (5…12) and calculates the rational approximation for all degrees: beginning from degree 2 to the highest possible degree. One repeats for every degree between 2 and the largest one the calculation of the Taylor expansion and the rational approximation. All poles and corresponding amplitudes for every degree are calculated and stored them in a master pole file. The latter is required for the stability analysis. Regarding this procedure, two comments are in order:

1. It is is very important to choose carefully the expansion point $`p_0`$ of the Taylor expansion. Once it is chosen the Taylor coefficients will be determined and from them the poles of the Laplace transform are obtained along with the corresponding amplitudes.
2. When the poles and amplitudes are found one has to assess the implications of the signs of the amplitudes and poles. For example: If the data shows a monotonic falloff, then the poles should all be real negative. If the data show oscillatory behaviour, then one must expect imaginary or complex poles. The poles and amplitudes stored in pole. file for all degrees calculated, from 2 to the highest one (5..12) will be analysed by a stability check procedure that will establish the history of occurrence of each pole as the degree has been varied.

### II.4 Pole stability analysis:

After poles and amplitudes are calculated from the Laplace and rational approximation, one has to study the number of occurrence of each pole with its corresponding amplitude. Since we are looking at the frequency of occurrence of poles in the two-dimensional complex plane, we have to look at small areas where a given pole keeps showing up persistently beyond a certain degree. For a given rough estimate of the size of the area in which the occurrence of poles is confined (an interactive test is provided for the size of the radius of the circular area). The sensitivity analysis of the presumed stable poles and their occurrence in that small area is performed in order to evaluate the robustness of the poles found.

### II.5 Consistency analysis:

As a final analysis of the soundness of the results, the software provides three capabilities:

1. The sequence of ratios of the Taylor coefficients should converge towards the smallest modulus of the poles found. This is a test of the accuracy of the calculation of the Taylor coefficients.
2. Once the stable amplitudes and poles are found and have passed all tests descibed above, one can define an analytical function from them that can generate data to be compared with the initial input data.
3. If the analysis fails or one wishes to perform a robustness analysis, it is required that the Taylor coefficients are well calculated since all the subsequent analysis is based on them. The software allows the selection of several integration techniques along with several possible interpolation routes for the data to analyse. One may compare the values of the coefficients under the different assumptions and choose the satisfactory ones.

## III Detailed description of the processing steps:

1. Preprocessing:
Performing the smoothing of the data, estimating the noise contained in it and altering the sampling rate at will for accuracy considerations during the calculation of the Taylor coefficients.
2. Finding the optimal expansion point $`p_0`$:
Two procedures are used:

1. The Laplace transform of the data is performed in order to find the approximate locations of the poles. The optimal expansion point $`p_0`$ should be larger than the largest pole found.
    2. Determination of the optimal point: This performs the Taylor expansion of the Laplace transform of the data and estimates the damping coefficient required for rapidly varying data values. The Taylor coefficients are plotted as a function of their order for the given degree and set of points chosen.
3. Rational approximation to the Taylor expansion of the Laplace transform of the data:
This step is to perform the rational approximation and find the poles and amplitudes for all degrees starting from 2 up to the highest degree chosen.
4. Stability Analysis:
This step performs the stability analysis. The stable poles and amplitudes are found and stored.
5. Consistency and final checks:

1. Generating an analytical function from the stable poles and amplitudes found. The accuracy of the Taylor coefficients is being granted by the graph showing whether their ratios approch the minimal pole modulus or not as their order increases.
    2. Generating data points from the constructed analytical function and adding synthetic noise. One might even add Gaussian or Uniform noise in order to make it as similar as possible to the initial data.
    3. Accuracy testing: This consists in obtaining the Taylor coefficients with various methods depending on the data available and the accuracy needed. It is advised to see how the various ways of calculating the coefficients alters the final results (poles and amplitudes).

## IV Procedure summary:

The summary reviews every step and indicates how the required input data is processed and how the results obtained at every step described previously are checked:

1. Preprocessing: The user has to provide the data file name, the number of data points and how the file is organized (data only, index and data or x, y data). In the filtering step, an estimate of the noise contained in the data as well as a smoothing parameter should be provided. The spline procedure returns an estimate of the noise and if the smoothing parameter is well chosen the initial estimate of the noise and the returned noise figure should be roughly in agreement. Finally it is possible to choose a different sampling rate of the data once it has been smoothed and splined. This has implications on the accuracy of the Taylor coefficients and a test exists in the software to check this accuracy.
2. Optimal expansion point: Assuming a range of points about which the Taylor coefficients are calculated for a chosen degree, the software returns the graphs and values of these coefficients versus the degree.
3. Rational approximation: Providing an upper value of the degree for which the analysis is performed, the software returns all poles and amplitudes for all degrees between degree 2 and the largest degree chosen. It returns also a consistency graph showing the sequence of ratios of Taylor coefficients for increasing order versus the minimal modulus of the poles found.
4. Stability Analysis: For a given radius (within bounds) of the circular neighborhood about which the frequency of occurrence of a given pole is examined, this analysis is done. The software returns the history of occurrences of the poles and determines the stable ones.
5. Consistency Analysis: This is done automatically through the display of the graph of the ratios of the Taylor coefficients versus order with respect to the smallest modulus of the poles found with the Padé-Laplace rational approximation. A new set of data from the stable poles and amplitudes found can be generated and compared it to the initial input data. It is advised to examine the sensitivity of the Taylor coefficients with respect to the various ways of calculating them, depending on the type of data available.

## V Examples:

### V.1 Synthetic example: Degree 12 Lanczos

We start first with a synthetic example. The parameters are:

Amplitudes:
3.00000
0.400000
3.40000
6.20000

Corresponding exponential arguments (poles):
(-1.0000000000000, 0.)
(-3.0000000000000, 0.)
(-5.0000000000000, 0.)

Damping parameter: -4.16

1. Determination of the optimal expansion point and the damping coefficient:
the taylor expansion coefficients c(r) are listed according to their absolute values whereas r!c(r) are given in algebraic values:

$`p`$ = 1.000D-01,

c(1)= 2.6760968316948
c(2)= -0.92274649709482
c(3)= 0.46139353620293
c(4)= -0.31918555072345
c(5)= 0.26204150623143
c(6)= -0.22997288909213
c(7)= 0.20656813173070
c(8)= -0.18701514280649
c(9)= 0.16977029842334
c(10)= -0.15425931881287
c(11)= 0.14021104327295
c(12)= -0.12745666366475
c(13)= 0.11586714821751
c(14)= -0.10533295150486

$$\begin{array}{cc}\hfill & c(1)\hfill \\ \hfill & c(2)\hfill \\ \hfill & c(3)\hfill \\ \hfill & c(4)\hfill \\ \hfill & c(5)\hfill \\ \hfill & c(6)\hfill \\ \hfill & c(7)\hfill \\ \hfill & c(8)\hfill \\ \hfill & c(9)\hfill \\ \hfill & c(10)\hfill \\ \hfill & c(11)\hfill \\ \hfill & c(12)\hfill \\ \hfill & c(13)\hfill \\ & c(14)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & 1!c(1)\hfill \\ \hfill & 2!c(2)\hfill \\ \hfill & 3!c(3)\hfill \\ \hfill & 4!c(4)\hfill \\ \hfill & 5!c(5)\hfill \\ \hfill & 6!c(6)\hfill \\ \hfill & 7!c(7)\hfill \\ \hfill & 8!c(8)\hfill \\ \hfill & 9!c(9)\hfill \\ \hfill & 10!c(10)\hfill \\ \hfill & 11!c(11)\hfill \\ \hfill & 12!c(12)\hfill \\ \hfill & 13!c(13)\hfill \\ & 14!c(14)\hfill \end{array}$$

Optimal expansion point $`p_0`$=0.4
2. Laplace and rational approximation calculations for degrees= 2…12
3. Stability study of the poles versus radius:
radius=1
number of stable poles= 8
pole # 1 (-4.0787128997673, 0.) appeared 11 times
pole # 2 (-1.0821008803446, 0.) appeared 14 times
pole # 3 (-3.0000000476863, 0.) appeared 10 times
pole # 4 (3.8482758680100, 0.) appeared 3 times
pole # 5 (-8.6468139578141, 0.) appeared 2 times
pole # 6 (1.8636510844379, 0.) appeared 6 times
pole # 7 (-6.3269596366764D-02, -0.66701505012643) appeared 26 times
pole # 8 (-6.3269590975568D-02, 0.66701507038499) appeared 5 times
the number of occurrences of some pole is larger than the highest degree = 12

Going to a smaller radius=10-2

Number of stable poles= 10

pole # 1 (-4.9999999672223, 0.) appeared 10 times
pole # 2 (-3.0000000476863, 0.) appeared 10 times
pole # 3 (-0.99999996721772, 0.) appeared 10 times
pole # 4 (3.8482758680100, 0.) appeared 2 times
pole # 5 (-8.6468139578141, 0.) appeared 2 times
pole # 6 (1.8636510844379, 0.) appeared 2 times
pole # 7 (-6.3269596366764D-02, -0.66701505012643) appeared 2 times
pole # 8 (-6.3269590975568D-02, 0.66701507038499) appeared 2 times
pole # 9 (2.7049633928318, -1.6749074721005D-09) appeared 2 times
pole # 10 (0.40001532007036, 2.9758199321717D-03) appeared 2 times
4. Consistency check:
Occurrence of the poles versus degree (1 for occurrence and 0 for no-occurrence):

pole # 1: 0 1 1 1 1 1 1 1 1 1 1
pole # 2: 0 1 1 1 1 1 1 1 1 1 1
pole # 3: 0 1 1 1 1 1 1 1 1 1 1
pole # 4: 0 0 1 0 0 0 0 0 1 0 0
pole # 5: 0 0 0 1 0 0 0 1 0 0 0
pole # 6: 0 0 0 1 0 0 0 1 0 0 0
pole # 7: 0 0 0 0 1 0 1 0 0 0 0
pole # 8: 0 0 0 0 1 0 1 0 0 0 0
pole # 9: 0 0 0 0 1 0 1 0 0 0 0
pole # 10: 0 0 0 0 0 0 1 0 0 0 0

After reducing radius to 10-5

Number of stable poles= 3
pole # 1 (-4.9999999672223, 0.) appeared 10 times
pole # 2 (-3.0000000476863, 0.) appeared 10 times
pole # 3 (-0.99999996721772, 0.) appeared 10 times
pole # 1: 0 1 1 1 1 1 1 1 1 1 1
pole # 2: 0 1 1 1 1 1 1 1 1 1 1
pole # 3: 0 1 1 1 1 1 1 1 1 1 1
Conclusion of the analysis:
Rational approximation degree= 12
Expansion point: 0.400000, Damping parameter= -4.16000
5. Stable poles and their corresponding amplitudes:

amplitude # 1 (6.1999983164819, -9.6235532915790D-10)
pole # 1 (-5.0000004239243, -1.1584249121758D-11)
amplitude # 2 (3.4000016676147, 1.2316899810646D-08)
pole # 2 (-3.0000003747023, 9.9268993397018D-12)
amplitude # 4 (0.39999985296666, -1.8860442419695D-07)
pole # 4 (-0.99999997587661, 2.7211949994764D-10)
for comparison, the input amplitudes and arguments are:

3

0.400000
3.40000
6.20000

(-1.0000000000000, 0.)
(-3.0000000000000, 0.)
(-5.0000000000000, 0.)

### V.2 Decaying only exponentials

5

(200.000,0.)
(1000.000,0.)
(500.000,0.)
(800.000,0.)
(300.000,0.)

( -50.000000000000, 0.)
( -5.0000000000000, 0.)
( -0.50000000000000, 0.)
( -0.25000000000000, 0.)
( -0.15000000000000, 0.)

Damping parameter= -5.00

Expansion point $`p_0`$=1.2, damping parameter=-5.

amplitude # 1 ( 200.44842041705, -6.4297285954932D-08)
pole # 1 ( -50.199847855334, 1.7266006542437D-09)

amplitude # 2 ( 1000.06460142389, -5.7036189955559D-06)
pole # 2 ( -5.0002454942104, -1.3064326900177D-08)

amplitude # 3 ( 499.40738201272, 1.3497166341200D-05)
pole # 3 ( -0.50015131820020, -1.6254860197388D-11)

amplitude # 4 ( 798.94082988133, 2.3349965623979D-05)
pole # 4 ( -0.25024465114282, 0.)

amplitude # 5 ( 301.65857589825, 9.1799870303999D-06)
pole # 5 ( -0.15016920303216, 0.)

amplitude # 6 ( -6.8688938206640D-08, -3.5775972803453D-08)
pole # 6 ( 0.51069582321812, 0.17487380533828)

amplitude # 7 ( -6.8686582230124D-08, 3.5777421059967D-08)
pole # 7 ( 0.51069582321812, -0.17487380533828)

amplitude # 8 ( 3.9828988025503D-05, -1.2554952329213D-05)
pole # 8 ( 2.0290173507300, 1.2045224512252)

amplitude # 9 ( -2.3797048952317D-05, -2.7704251925653D-05)
pole # 9 ( 2.0290173521374, -1.2045224398712)

amplitude # 10 ( 5.4490417778712D-07, 3.4864743102636D-14)
pole # 10 ( 2.6347983285165, 0.)

Rational approximation degree= 10

Damping parameter= -5.00, Expansion point $`p_0`$: 1.20

Stable poles:

amplitude # 1 = ( 200.44842041705, -6.4297285954932D-08)
pole # 1 = ( -51.104140221450, 0.)

amplitude # 2 = ( 1000.06460142389, -5.7036189955559D-06)
pole # 2 = ( -5.3490716102094, 0.)

amplitude # 3 = ( 499.40738201272, 1.3497166341200D-05)
pole # 3 = ( -0.50758978101004, 0.)

amplitude # 4 = ( 798.94082988133, 2.3349965623979D-05)
pole # 4 = ( -0.26205674432617, 0.)

amplitude # 5 = ( 301.65857589825, 9.1799870303999D-06)
pole # 5 = ( -0.15111811290477, 0.)

### V.3 Oscillating and decaying exponentials

5
(400.000,0.)
(175.000,0.)
(175.000,0.)
(-500.000,0.)
(250.000,0.)
( -2.0000000000000D-03, 0.)
( -8.0000000000000D-03, 3.0000000000000D-02)
( -8.0000000000000D-03, -3.0000000000000D-02)
( -2.0000000000000D-02, 0.)
( -0.20000000000000, 0.)

Stability radius = 10-2
Number of stable poles = 5

Stable poles and amplitudes:
amplitude # 1 = ( 250.00011450201, 9.2737811992658D-07)
pole # 1 = ( -0.20001227827931, 0.)

amplitude # 2 = ( -499.99976160965, -8.0435000145924D-04)
pole # 2 = ( -2.0047256715262D-02, 4.3624394471021D-13)

amplitude # 3 = ( 175.00010854502, -8.5323335360501D-05)
pole # 3 = ( -8.0240426328942D-03, -2.9980187019490D-02)

amplitude # 4 = ( 175.00006247227, 2.3514996070493D-04)
pole # 4 = ( -8.0240426329283D-03, 2.9980187019495D-02)

amplitude # 5 = ( 399.99945705790, 6.6704669352807D-04)
pole # 5 = ( -1.9888744154064D-03, 0.)

### V.4 Experimental example borrowed from fast spectroscopy

The example below Fig.1 shows an experimental noisy signal borrowed from fast laser spectroscopy. After preprocessing the data and estimating the value of the noise present in the data, we perform the Padé-Laplace analysis and the results are displayed in Fig.2. The difficulty in fitting this example is due not only to the presence of noise but also of the large initial value in the data. Despite these difficulties, G-MAPLE was able to extract the exponential arguments and corresponding amplitudes successfully for a rough value of $`p_0`$. Slightly changing $`p_0`$ gave a minor improvement of the fitting. The obstacle in finding a better fit stems from the fact it is harder to find a good expansion point $`p_0`$ for real noisy data. The challenge is to find a better automatic procedure for the determination of $`p_0`$ in the case of noisy data. In that case, least-squares cadzow  or fast least-squares estimation methods provide better answers demeure .

Acknowledgements
The author wishes to acknowledge helpful discussions with Randy Borle and Art Monk that helped improve substantially the graphical user interface of the software. He also thanks Daniel Houde for kindly providing two fast spectroscopy data sets for testing the software, and P. Rasolofosaon for providing reprints/preprints of his work.


text:
# On emission-line spectra obtained from evolutionary synthesis models

## 1 Introduction and motivation

In recent years a great deal of effort has been put into modeling more realistic ionising fluxes of massive stars (e.g., the wm-basic code by Pauldrach et al. 2001, and references therein). These new model atmospheres, as included in evolutionary synthesis codes, are expected to produce more reliable ionising spectra, which in turn should allow a better interpretation of the emission line spectra of ionised nebulae, through the use of photoionisation codes (e.g., Smith et al. 2002).

Yet the ionising spectra of stellar clusters are produced by the most massive stars, whose absolute numbers can be strongly affected by the incomplete sampling of the underlying Initial Mass Function (IMF). The effects of incomplete sampling on the observed properties of star-forming regions can be currently evaluated in two main ways: either by means of Monte Carlo synthesis models, or by means of an appropriate statistical formalism applied to the results of analytical synthesis models. By analytical synthesis models, or analytical models for short, we will in the following refer to those population synthesis models computed by codes that neglect the issue of the incomplete sampling, and assume instead that the mass spectrum of stellar clusters can be represented by a continuous, analytical function that coincides with the underlying IMF. On the other hand, Monte Carlo codes simulate stellar populations by the stochastic generation of stars, which is stopped when either a given number of stars or a given cluster mass is reached; the mass probability function for each star is given by the IMF, in such a way that the actual mass spectrum approaches the IMF as the number of stars increases, whereas large deviations occur for small clusters. This is assumed to be a more realistic representation of real stellar clusters than the one provided by analytical codes.

The objective of this work is to study quantitatively the impact on emission line spectra of the sampling effects in the IMF, and to evaluate its consequences in the determination of the physical properties of H ii regions and galaxies. This is not a trivial problem. A simple, perhaps naïve, approach would go through the following steps: (i) building synthetic clusters with different numbers of stars with a stochastic sampling of the IMF, ensuring in this way that only integer numbers of stars are included; (ii) obtaining the corresponding Spectral Energy Distributions (SEDs); (iii) using the SEDs as inputs to a photoionisation code; and (iv) drawing conclusions based on the overall results. However, since the number of free parameters is pretty large, it is far more useful to split the problem into several smaller steps, each of which addressing a particular source of uncertainty. A possible procedure adopting this strategy is the following:

* A first step is to evaluate quantitatively the sampling effects on the ionising spectrum. Although the computation of photoionisation models is not involved in this step, it is important to characterize the resulting ionising flux with a view on the physics of photoionisation.
* A second step is to determine how the different emission lines scale with the size (mass/number of stars) of the stellar cluster. For example, the intensity of hydrogen recombination lines roughly scale with the number of ionising photons above the threshold energy, and hence with the number of massive stars in the cluster (for a given cluster age and given properties of the emitting gas). However, a similar scaling relation -provided it exists- has not yet been established for the case of forbidden emission lines.
* A third step is to investigate the influence of both the stellar and the nebular properties on the results of photoionisation models. A large grid of photoionisation models would be required to explore exhaustively the parameter space, both along stellar population dimensions (i.e., IMF slope and mass limits, ages, stellar metallicity, model atmospheres, etc.), and along gas properties dimensions (density profile, chemical abundances, radius of the nebulae, covering factor, etc.). In this step we can partially rely on extensive studies performed by different groups (e.g. Dopita et al. 2000, and references therein). Feedback from observations must also be taken into account.

In this series of papers, we address each of these problems in turn. This first paper presents a quantitative evaluation of sampling effects on the predicted ionising flux for a wide range of input parameters in synthesis models, allowing us to obtain some conditions necessary, but not sufficient, to explain the correlations of observational data in emission-line diagnostic diagrams. In addition, we determine a mass limit under which stellar clusters should be better modeled with single individual stars rather than with stellar clusters. In a second paper we will determine the scaling relations of the emission-line spectrum with the initial mass/number of stars in the cluster. The third, and last, paper of this series will be devoted to the analysis of a set of Monte Carlo simulations of evolutionary synthesis models linked with a photoionisation code, where some specific properties will be assumed for the emitting gas. Throughout all these papers we will discuss the relationships between sampling effects and the ionising continuum, along with the resulting emission-line spectra.

The importance of observational feedback in this analysis cannot be overstated. A long-standing puzzle is posed by the positions in emission-line diagnostic diagrams of observed H ii regions ionised both by single stars and by stellar clusters cover an area narrower than the one expected on the basis of extensive grids of photoionisation models. For example, Dopita et al. (2000) built a grid of zero-age stellar clusters varying the metallicity and the ionisation parameter of the model nebulae, and found that the area covered on diagnostic diagrams by these models is much broader than the area covered by observational data, suggesting the existence of a hidden parameter that constrains the possible positions of real H ii regions in diagnostic diagrams. This problem has also been thoroughly discussed by Stasińska & Izotov (2003); these authors, who also consider the evolution of the ionising flux due to the aging of the stellar population, effectively highlight the difficulties of reproducing simultaneously all the observational constraints by means of traditional photoionisation models.

Since the essence of this problem is finding a way of reducing the predicted dispersion, one would naïvely expect that the problem would worsen when the IMF sampling effects are taken into account: that is, that the positions occupied in diagnostic diagrams by model nebulae would show a larger dispersion as a consequence of taking into account the dispersion in the ionising flux predicted by stochastic models. Instead, we show here that the intensities of the predicted flux at different energies are strongly correlated, a fact implying that no additional dispersion will result from the incomplete sampling of the IMF.

The structure of this paper is the following: in Section 2 we review briefly the evaluation of sampling effects and the definition of the Lowest Luminosity Limit. The application of these concepts to the case of ionising spectra is presented in Section 3. Our main results are explained in Section 4, and the conclusions are presented in Section 5.

## 2 Consequences of sampling effects in synthesis models

In this section we describe how sampling effects in the IMF are evaluated. Some of these results have either been presented in previous papers, or can be found in books on statistics, e.g. Kendall & Stuart (1997): nevertheless, we have chosen to summarize them here rather than in an appendix, because they are essential for the forthcoming discussion. Readers that are already familiar with these tools may jump directly to Section 3.

### 2.1 The Lowest Luminosity Limit

Sampling effects produce a dispersion in the results of synthesis models, and should therefore be taken into account when the properties of real clusters are interpreted by means of synthesis models. The relevance of such effects depends on the size of the system studied. A starting point to traduce this statement into quantitative terms is the definition of a criterion, related to the cluster size, allowing one to establish whether sampling effects in a particular cluster are playing a fundamental role or not. This alternative in turn would condition the mode of application of synthesis models: clusters small according to this criterion should necessarily be modeled with full consideration of sampling effects, whereas larger clusters could be modeled by means of traditional analytical methods.

A suitable criterion to classify systems in this way is provided by the Lowest Luminosity Limit (or LLL: Cerviño & Luridiana 2003a, b). The LLL is defined as the luminosity of the most luminous individual star compatible with the physical properties (age, metallicity, etc.) of the cluster, and the criterion to follow is that an observed cluster whose luminosity is below the LLL cannot be modeled by means of a method that neglects sampling effects. The logic underlying this definition is that an analytical model computed to reproduce a cluster below the LLL necessarily includes fractions smaller than one of the relevant stars, and therefore lack physical meaning. Note that the restriction imposed by this criterion is quite loose, since it constrains only those situations in which strong statistical effects are undoubtedly present; but important statistical effects may well be present even in clusters above the LLL.

Since the average luminosity $`A`$ of a cluster scales linearly with the total mass $``$, the LLL is naturally associated to a lower mass limit $`_A^{min}`$, which is the mass of a stellar system with a completely sampled IMF and luminosity equal to the LLL.

The numerical value of the LLL depends on the isochrones and the atmosphere models used, does not depend on the IMF, and is only weakly dependent on the star formation history; on the contrary, it depends strongly on the age and metallicity of the synthetic cluster. Finally, the LLL obviously depends on the energy band in which the luminosity is defined. For a given luminosity band $`A`$ the LLL defines an intrinsic limit below which a computation of the statistical dispersion associated to the incomplete sampling of the IMF is unavoidable for a meaningful application of synthesis models. Above the LLL, it is possible to provide a rough estimation of the dispersion and a guideline to determine the relevance of sampling effects: for example, Cerviño & Luridiana (2003a) estimate that the dispersion in the results of synthesis models for clusters with initial masses about $`10\times _A^{min}`$ is equal or larger than 10% in the optical and infrared bands. Below the LLL, a proper statistical formalism is required to obtain quantitative estimations of the relative dispersion in the predicted quantities; this topic will be covered in the next section. However, note that below the LLL the predicted average values of observables that do not scale linearly with $``$, like logarithmic or rational functions of luminosities (i.e., equivalent widths or colours) can be severely biased with respect to actual observations (see Cerviño & Valls-Gabaud 2003), so that even a sophisticated statistical formalism can fail to produce meaningful predictions.

### 2.2 Estimation of sampling effects

As stated in Section 1, sampling effects in synthesis models can be directly estimated by Monte Carlo methods, or alternatively by a statistical formalism applied to the results of analytical synthesis models. A formalism of this kind, based on the original formalism proposed by Buzzoni (1989), has been developed in recent years by Cerviño et al. (2000, 2001, 2002b); and Cerviño & Valls-Gabaud (2003). The method can be briefly summarised as follows. Let us assume that $`w_i`$ is the number of stars within a given mass range, normalised to the total mass; this value is given by the IMF and the star formation history. Assuming a Simple Stellar Population model (SSP, i.e. a model with coeval star formation, also called Instantaneous Burst model), $`w_i`$ is given by integrating the IMF over the corresponding mass range. Each $`w_i`$ is approximately the mean value of a Poisson distribution, and so the variance $`\sigma _i^2`$ of each $`w_i`$ is equal to $`w_i`$. This assumption is a good approximation if the mass interval that defines the $`w_i`$ values is narrow enough: see Cerviño et al. (2002b) for a complete discussion of this topic.

Given an observable property $`a_i`$ of a given star $`i`$, its contribution to the average normalised integrated property $`\mu _A^{ssp}`$ obtained by the synthesis model is given by $`w_ia_i`$, with a variance $`\sigma _i^2a_i^2=w_ia_i^2`$. The total variance of the average observable $`\mu _A^{ssp}`$ is the sum of all the variances. Hence the relative dispersion is simply

$$\frac{\sigma _A^{ssp}}{\mu _A^{ssp}}=\frac{(_iw_ia_i^2)^{1/2}}{_iw_ia_i}=\frac{1}{\sqrt{𝒩_A^{ssp}}},$$
(1)
where the last term gives the definition of $`𝒩_A`$ introduced by Buzzoni (1989). Since $`w_i`$ scales with the initial mass $``$ transformed into stars in the cluster, the effective number of stars contributing to observable $`A`$ also scales with the initial mass. This implies that, if the $``$ value of an observed cluster is known, its predicted average $`A`$ value is $`\mu _A=\mu _A^{ssp}\times `$, and the dispersion around $`\mu _A`$ is given by $`𝒩_A=𝒩_A^{ssp}\times `$.

In addition, given another property $`b_i`$, with contribution to the integrated average property $`\mu _B^{ssp}`$ given by $`w_ib_i`$ and variance $`(\sigma _B^2)^{ssp}`$, it is also possible to obtain the linear correlation coefficient between the observables $`A`$ and $`B`$:

$`\rho (A,B)`$ $`=`$ $`{\displaystyle \frac{\mathrm{cov}(A,B)}{\sigma _A\sigma _B}}={\displaystyle \frac{\mathrm{cov}(A,B)^{ssp}}{\sigma _A^{ssp}\sigma _B^{ssp}}}`$ (2)
$`=`$ $`{\displaystyle \frac{_iw_ia_ib_i}{\sqrt{(_iw_ia_i^2)(_iw_ib_i^2)}}},`$
where $`\mathrm{cov}(A,B)`$ is the covariance. Although $`\mathrm{cov}(A,B)`$ scales with $``$, $`\rho (A,B)`$ is independent of it. The correlation coefficient can also be obtained by a linear regression of $`A`$ vs. $`B`$: a value of $`\rho (A,B)=1`$ implies that the two quantities are linearly dependent. In terms of a hypothetical diagnostic diagram of $`A`$ vs. $`B`$, $`\rho (A,B)1`$ implies that the data points follow a narrow sequence.

## 3 Methodology

The simplest way to characterize the ionising flux of a stellar population is to consider both its intensity (in number of photons emitted) and its shape (in terms of “effective rate of ionising photons” defined below). The emission rate of photons capable of ionising a particular ion ($`Q(\mathrm{H}^0)`$, $`Q(\mathrm{He}^+)`$, etc.) with ionisation threshold frequency $`\nu _0`$ is given by:

$$Q(X^i)=_{\nu _0}^{\mathrm{}}\frac{L_\nu }{h\nu }𝑑\nu .$$
(3)
As for the shape of the ionising flux, it is important to note that not all ionising photons are equivalent since the photoionisation cross section depends strongly on $`\nu `$ (see Fig. 1). The effectiveness of photons in ionising a ion $`X^i`$ is measured by the balance between $`X^i`$ and $`X^{i+1}`$ at a given point in the nebula:

$$\frac{N(X^{i+1})N_e}{N(X^i)}=\frac{_{\nu _i}^{\mathrm{}}(4\pi J_\nu /h\nu )a_\nu (X^i)𝑑\nu }{\alpha _G(X^i,T)},$$
(4)
where $`N(X^{i+1})`$, $`N(X^i)`$, and $`N_e`$ are the number densities of ions $`X^{i+1}`$, $`X^i`$, and electrons respectively, $`\alpha _G(X^i,T)`$ is the recombination coefficient of the ground level of $`X^{i+1}`$ to all levels of $`X^i`$, $`a_\nu (X^i)`$ is the absorption cross section of the ion $`X^i`$ from its ground level with threshold $`\nu _i`$, and $`4\pi J_\nu `$ is the local flux density (Osterbrock 1989). Since we are ultimately interested in determining the emission spectra, which in turn depend on the relative ionic abundances, it seems physically meaningful to weigh the emission rate of photons capable to ionise a particular ion by the photoionisation cross section for that ion, in order to describe the luminosity actually seen by each ion; The resulting quantity, which we call the “effective rate of ionising photons” is defined by the expression:

$$Q_{a_\nu }(X^i)=_{\nu _i}^{\mathrm{}}\frac{L_\nu }{h\nu }a_\nu (X^i)𝑑\nu .$$
(5)
This quantity is obviously not an observable, but it will help to perform first-order estimations of the sampling effects on the ionising flux, via Eq. 4. We have computed $`Q_{a_\nu }(X^i)`$ for H0, He0, He+, N0, N+, N++, O0, O+, O++, S0, S+, S++, C0, Ne0, Ne+, Ar0, Ar+, and Ar++, although in this work the Ar, Ne, and C ions will not be discussed111The complete set of results can be found at the CMHK web pages at http://laeff.inta.es/users/mcs/SED/.. The corresponding $`a_\nu `$ values have been obtained from the subroutine phfit.f used for this purpose in the photoionisation code Cloudy (Ferland 1996), and are shown in Fig. 1. It is important to recall that since the recombination lines of $`X^i`$ arise from the recombination of $`X^{i+1}`$, they depend on $`Q_{a_\nu }(X^i)`$ (e.g., $`Q_{a_\nu }(\mathrm{H}^0)`$ is related to the H i lines). On the other hand, since forbidden lines of $`X^i`$ arise from $`X^i`$ levels, they depend on $`Q_{a_\nu }(X^{i1})`$ (e.g. $`Q_{a_\nu }(\mathrm{O}^+)`$ is related to the \[O iii\] lines, and not to the \[O ii\] lines).

## 4 Results

In this section we will describe our results. All the computations assume an instantaneous burst of star formation and a Salpeter (1955) IMF in the mass range 0.1 to 120 M.

To assess the possible systematic effects associated with the choice of stellar tracks, model atmospheres, interpolation schemes, etc, we have used two different evolutionary synthesis codes:

* The code by Cerviño et al. (2002a) (hereinafter CMHK), which adopts evolutionary tracks with standard mass-loss rates (Schaller et al. 1992; Schaerer et al. 1993a, b; Charbonnel et al. 1993) and the model atmospheres by Schaerer & de Koter (1997) (CoStar) for main sequence hot stars more massive than 20 M, by Schmutz et al. (1992) for WR stars, and by Kurucz (1991) for the remaining stars.
* A modified version of Starburst99 (Leitherer et al. 1999), where the main modification is the computation of the $`Q_{a_\nu }(X^i)`$ values and the corresponding minimum masses. The use of Starburst99 is only for comparison purposes, since we are interested in studying the impact of the new atmosphere models by Smith et al. (2002) for WR and hot stars, implemented in this code. Note that Starburst99 assumes Blackbody spectra for stars with $`\mathrm{log}T_{eff}>4.778`$ or $`\mathrm{log}T_{eff}<3.3`$ except for WR stars (defined by the stellar temperature, $`\mathrm{log}T_{}>4.4`$, the hydrogen surface abundance $`X_s<0.4`$, and the minimum mass for WR formation). The stellar tracks used have high mass-loss rates (Meynet et al. 1994).

In both codes, the spectrum associated to a given star is given by the closest model atmosphere found either in the $`\mathrm{log}gT_{eff}`$ plane or in the $`gT_{eff}`$. For the computation of the isochrones, both codes follow the prescriptions described in Cerviño et al. (2001).

### 4.1 Effective rate of ionising photons $`Q_{a_\nu }(X^i)`$

The time evolution of $`Q_{a_\nu }(X^i)`$ predicted by Starburst99 and the CMHK code is shown in Fig. 2 and 3 respectively, where the five main panels correspond to the metallicities indicated222Note the different scale along the $`Q_{a_\nu }`$-axis for $`Z`$=0.020 and $`Z`$=0.040.. Each main panel is further subdivided into four sub-panels, each showing the results for a different set of ions, coded according to the lines and symbols indicated. The top-left panel corresponds to H0, He0, and H+; the top-right panel to N0, N+, and N++; the bottom-left panel to O0, O+, and O++; and the bottom-right panel to S0, S+, and S++.

There are striking differences between the $`Q_{a_\nu }(X^i)`$ behaviour predicted by the two codes, a fact that emphasises the importance of the models’ input – evolutionary tracks and atmosphere models – as a source of systematic effects. An example is the difference in the bumps of the $`Q_{a_\nu }(\mathrm{He}^+)`$ curves: these bumps reveal the life cycle of WR stars, which provide a substantial amount of hard photons. The low $`Z`$ models computed with Starburst99 show such bumps at about 3 Myr, whereas the CMHK models do not; the difference is mainly due to the different evolutionary tracks, because enhanced mass-loss rates produce larger numbers of WR stars. A further difference is the larger $`Q_{a_\nu }(\mathrm{He}^+)`$ value at high $`Z`$ in the CMHK code as compared to Starburst99. This difference can be ascribed to differences in the atmosphere models implemented in the two codes: in general, the CMHK code produces a harder ionising flux due to the use of the CoStar and Schmutz et al. (1992) atmosphere models. This effect is larger for high metallicities during the WR phase. The impact of different atmosphere models in synthesis codes has been recently discussed by Smith et al. (2002), and we refer to that paper for further details.

There are two other interesting features in Fig. 2 and 3 that deserve discussion:

1. Although the behaviours of $`Q_{a_\nu }(\mathrm{He}^+)`$ and $`Q_{a_\nu }(\mathrm{O}^{++})`$ are remarkably similar, as a consequence of their ionisation edges being almost identical (54.418 eV for He+ and 54.936 eV for O++), there are slight differences between them due to the greater sensitivity to low energy photons of $`a_\nu (\mathrm{He}^+)`$ with respect to $`a_\nu (\mathrm{O}^{++})`$. It follows that $`Q_{a_\nu }(\mathrm{O}^{++})`$ decreases more rapidly than $`Q_{a_\nu }(\mathrm{He}^+)`$ when the hard ionising flux decreases, that is for ages $`t>4`$ Myr.
2. The tails of the $`Q_{a_\nu }(\mathrm{He}^+)`$ and $`Q_{a_\nu }(\mathrm{O}^{++})`$ curves during the post-WR phases show an interesting behaviour in both sets of results. For example, at $`Z`$=0.008 there is a $`Q_{a_\nu }(\mathrm{He}^+)`$ peak around 6 Myr in the case of stellar tracks with enhanced mass-loss rates (Fig. 2), and a peak around 6.5 Myr in the case of stellar tracks with standard mass-loss rates (Fig. 3). Yet, although the bumps in these curves reveal the presence of WR stars in the stellar population, the WR phase ends at 5.9 Myr in the case of enhanced mass-loss rates, and at 4.7 Myr in the case of standard mass-loss rates. The explanation of this apparent contradiction lies in the track interpolation technique used to obtain the isochrones: even though it is assumed that the evolutionary tracks follow a continuous sequence, there is a discontinuity in the stellar tracks between WR and non-WR phases. The situation is illustrated in Figure 4, where isochrones for 4.5, 6 and 6.3 Myr obtained from Starburst99 are compared with the evolutionary tracks for 60, 40 and 25 M with enhanced mass-loss rates. Note that some stars in the 6 Myr isochrone (where no WR stars should be present) reach $`\mathrm{log}T_{eff}`$ larger than 4.4, and they are therefore (wrongly) assigned to the WR type. A similar problems occurs in the CMHK code, and in general in all synthesis codes, since continuity in the tracks is implicitly assumed (see, however, an alternative approach presented in Cerviño & Mas-Hesse 1994). The isochrones are mathematically consistent, but this is not a guarantee that they are physically consistent. This problem is amplified in the tracks with enhanced mass-loss rates since adjacent tracks are far more different than the corresponding tracks with standard mass-loss rates. This may superficially be considered a technical detail (Maeder 2002), but since it may produce unphysical results, it is worth investigating it in detail (Cerviño 2003, in preparation).

### 4.2 Minimum masses $`^{min}[Q_{a_\nu }(X^i)]`$

In analogy to our definition of the minimum mass associated to the LLL, we will define in the following a minimum mass associated to the effective rate of ionising photons, $`^{min}[Q_{a_\nu }(X^i)]`$, which is defined as the mass of a stellar system with a completely sampled IMF and total effective rate of ionising photons equal to the maximum of all the effective rates provided by the individual stars in the system.

The evolution of $`^{min}[Q_{a_\nu }(X^i)]`$ for the ions selected before, computed with Starburst99 and the CMHK code, are shown in Fig. 5 and 6 respectively. The organisation of the figures is the same as in Fig. 2. Perhaps the most striking feature of these results is the huge mass range covered, which spans nearly three orders of magnitude, from about $`3\times 10^3`$ to more than 106 M for some of the important ions, independently of the code used. As expected, the larger the ionising potential, the larger the minimum mass $`^{min}[Q_{a_\nu }(X^i)]`$, since incomplete sampling affect mainly the most massive stars, which are those contribute most to this property. Likewise, S0 has the lowest $`^{min}[Q_{a_\nu }(X^i)]`$ value: the ionising potential of S0 is only 10.36 eV, so the effective rate of ionising photons for this ion is scarcely affected by sampling effects. On the other hand, He+ is the ion more affected by sampling effects: this result confirms that He ii lines cannot be used as a strong constrain on the properties of stellar clusters, as suggested by Luridiana et al. (2003b).

The most prominent peak in $`^{min}[Q_{a_\nu }(X^i)]`$ corresponds to the WR phase of the cluster. In the case of Starburst99, there is a tendency to obtain larger $`^{min}[Q_{a_\nu }(X^i)]`$ values during the WR-phase at intermediate metallicities ($`Z`$=0.008), and an asymmetric “U-shaped” curve at low metallicities. We have tested this behaviour obtaining the corresponding $`^{min}[Q_{a_\nu }(X^i)]`$ with the CMHK code using the evolutionary tracks by Meynet et al. (1994)333Which results are also available at the CMHK web pages., and the “U-shaped” behaviour remains. However, $`^{min}[Q_{a_\nu }(X^i)]`$ tends to be lower for larger metallicities. For the case of the CMHK code, the effect is explained easily: since the WR atmosphere models are the same for all the metallicities, the larger the metallicity, the larger the mass-loss rate, the lower the luminosity, and the lower $`^{min}[Q_{a_\nu }(X^i)]`$. In the case of Starburst99 the effect of the metallicity-dependent WR atmosphere models must also be taken into account together with the variation of the luminosity with metallicity. In fact, examining the ionising fluxes quoted for the WR stars in tables 3 and 4 of Smith et al. (2002), it is found that the flux of ionising photons is not a monotonic function of the metallicity: the largest values for the rate of ionising photons correspond to some of the WR models at $`Z`$=0.008 ($`Z`$=0.4 $`Z_{}`$ in their paper). This would be the origin of the larger $`^{min}[Q_{a_\nu }(X^i)]`$ value at this metallicity. This also shows the impact of the use of metallicity-dependent atmosphere models for WR stars.

The “U-shaped” behaviour during the WR-phase is due to the evolutionary tracks used: since the lower the number of stars that dominate the luminosity, the larger the $`^{min}[Q_{a_\nu }(X^i)]`$ value, its value depends on the relative rate of ionising photons emitted by WR stars with respect to the other hot stars in the synthetic cluster. In the tracks by Meynet et al. (1994), the WR stars at the beginning of the WR phase are also the most luminous stars in the cluster (see, for example, the 60 M track in Fig. 4) with the harder ionising flux: this fact produces the first peak in $`^{min}[Q_{a_\nu }(X^i)]`$. When the system evolves, the luminosity of WR stars decreases and the other non-WR, hot stars give only a small contribution (see, for example, the 40 M track in Fig. 4), leading to the decrease in $`^{min}[Q_{a_\nu }(X^i)]`$. However, since the non-WR stars become cooler as the age increases, at the end of the WR-phase WR stars become again the almost only source of ionising flux, producing the secondary peak. Note also that, for $`Z`$=0.040, the end of the WR-phase is marked by the highest peak in the curve: this behaviour is due to the longer duration of the WR-phase at this metallicity and to the dependency of the lifetime of massive stars with mass, which in these tracks is not monotonic: stars with initial mass of 85 and 120 M stars end their evolution as WR at 5 and 8 Myrs respectively. As a consequence, in such an age range and at this metallicity, the ionising flux is dominated by the less populated tail of the IMF, so that sampling effects dominate and a larger $`^{min}[Q_{a_\nu }(X^i)]`$ value is obtained.

The same general explanations also apply to the case of the CMHK results with standard mass-loss rate. However, for these tracks WR and non-WR stars have more similar luminosities and the first peak does not appear (the curve is “J-shaped”). In this case the maximum in $`^{min}[Q_{a_\nu }(X^i)]`$ occurs at $`Z`$=0.004, because the WR phase is very short and almost absent at $`Z`$=0.001. In fact, at this metallicity, there no star ends its evolution as a WR; the WR-phase only appears in the middle of the 120 M track, which ends as a supergiant.

There are also several peaks in $`^{min}[Q_{a_\nu }(X^i)]`$ at the oldest ages considered in these plots. These features are an artefact of the way in which model atmospheres are assigned to stars in the synthetic HR diagram: since the smooth evolutionary path of each star in the HR diagram is coupled to a discrete ensemble of atmosphere models, a peak is produced when stars suddenly switch from one atmosphere model to the next. In the case of Starburst99, the effect is not too strong because this code uses a finer atmosphere grid. In the case of the CMHK code, the switch from the CoStar to the Kurucz (1991) atmosphere models for stars in the main sequence around $`M=20`$ M produces the additional bumps in $`^{min}[Q_{a_\nu }(X^i)]`$ with the most prominent ones ending at $`t`$ 9.4 Myr for $`Z`$=0.001 and $`Z`$=0.004, $`t`$ 9 Myr for $`Z`$=0.008, $`t`$ 8.1 Myr for $`Z`$=0.020, and $`t`$ 7 Myr for $`Z`$=0.040, corresponding to the turn-off age of a 20 M star (point 13 in the tracks), the phase at which CoStar models are no longer used. Additional small bumps are found at the turn-off age of a 25 M. We want to remark again that these bumps are model artifacts, and that they would not appear if the atmosphere models were interpolated, instead of being computed by means of the “closest atmosphere model” approach.

Note that $`^{min}[Q_{a_\nu }(X^i)]`$ corresponds to the amount of gas transformed into stars, and not to the total mass of the cluster (stars plus gas, that is dynamical mass, which would be approximately a factor of ten larger), and that a Salpeter IMF in the mass range 0.1 to 120 M has been assumed. Table 1 gives conversion factors for different mass limits and also the numbers of stars for a Salpeter IMF. The results of Starburst99 are usually normalised to $`10^6`$ M in the mass range 1–120 M, while CMHK normalise the results to 1 M in the mass range 2–120 M.

Using the values in Table 1, a value of $`^{min}[Q_{a_\nu }(X^i)]`$ around $`3\times 10^3`$ implies that there are, on average, two stars in the mass range 40–120 M, or 21 stars in the mass range 8–120 M. Since the most massive stars are also the most luminous ones, the ionisation of the cluster is accounted for by a few stars. The presence of a few WR stars, with a harder ionising flux, produces the increases in $`^{min}[Q_{a_\nu }(X^i)]`$, especially for ions with higher energy edges. A straightforward conclusion follows from this result: the ionisation structure for clusters with $``$ close to $`^{min}[Q_{a_\nu }(X^i)]`$ would be better reproduced by the ionising spectrum of a single star rather than by the synthetic spectrum produced by a synthesis code, which always includes an ensemble of stars.

It is necessary to recall, in this context, that $`^{min}[Q_{a_\nu }(X^i)]`$ is, in fact, a rather restrictive lower limit for the use of integrated ionising spectra obtained by synthesis models. A better way to estimate when a composite spectrum begins to be a better approximation is given by the relative dispersion in $`Q_{a_\nu }(X^i)`$ corresponding to $`=^{min}[Q_{a_\nu }(X^i)]`$, which is shown in Fig. 7. As expected, the dispersion varies strongly with the age, metallicity and ion under consideration, but, in general, it has a value around 60% (or, equivalently, $`𝒩3`$). This means that in order to obtain a dispersion smaller than 10% ($`𝒩=100`$), clusters with initial masses $`>33\times ^{min}[Q_{a_\nu }(X^i)]`$ are needed. This is a good lower limit to ensure that the use of the composite ionising spectrum obtained from synthesis models is appropriate. For clusters with masses in the range \[$`^{min}`$, $`33\times ^{min}`$\] the situation is fuzzier, and any result obtained from codes which use a fully-sampled IMF must be taken with extreme caution. The only proper alternative solution in this case is to use Monte Carlo simulations.

### 4.3 Correlations between ionising bands

In this section, we will discuss the information that can be obtained through the study of the correlation coefficients $`\rho [Q_{a_\nu }(X^i),Q_{a_\nu }(Y^j)]`$ between the effective rate of photons of the ions $`X^i`$ and $`Y^j`$. For the sake of simplicity, we will use in the following the notation $`\rho (X^i,Y^j)\rho (Q_{a_\nu }(X^i),Q_{a_\nu }(Y^j))`$.

As indicated in Section 2.2, a linear correlation coefficient $`\rho (X^i,Y^j)`$ close to unity indicates that $`Q_{a_\nu }(X^i)`$ and $`Q_{a_\nu }(Y^j)`$ are essentially produced by the same stars within a stellar population. Note that the correlation between observable quantities is not affected by sampling effects in the IMF; nevertheless, we take advantage of the statistical formalism introduced in Section 2.2 for the analysis of the effects of incomplete sampling to compute the values of the correlation coefficients.

The correlation coefficients between effective rates of ionising photons may provide useful insights into the concept of ionisation correction factors (icfs) used in abundance determinations. A correlation coefficient close to one is a necessary condition to relate the abundance of two different ions in a ionised nebula independently of the ionising stellar population. The assumption of the existence of general recipes relating the abundances of two ions is conceptually analogous to the one underlying the definition of the icfs. However, a correlation coefficient close to one is not a sufficient condition for a correct definition of the icfs, since the properties of the gas might also be relevant in the determination of the ionisation structure of a nebula (Luridiana et al. 2003a).

The correlation coefficients between the effective rates of ionising photons of different ions may also be used to make predictions on diagnostic diagrams. However, in the case of diagnostic diagrams an additional difficulty appears, because they do not relate absolute intensities but ratios. Hence, e.g., the correlation coefficients $`\rho (\mathrm{H}^0,\mathrm{N}^0)`$ and $`\rho (\mathrm{H}^0,\mathrm{O}^+)`$ are not sufficient to derive a correlation in the \[N ii\] 6584/H$`\alpha `$ vs. \[O iii\] 5007/H$`\beta `$ diagram, since it is also necessary to know $`\rho (\mathrm{N}^0,\mathrm{O}^+)`$. In addition, even in the case of a full linear correlation ($`\rho =1`$), the tilt of the regression line and the position of the mean value are required, so a photoionisation model is always needed. For the time being, we will only obtain the correlation coefficients, to see whether correlations exist. In Paper III we will apply these results to make predictions on the dispersion in the diagnostic diagram.

The results of the calculation of the correlation coefficients are shown in Fig. 8 for selected ions and diagnostic diagrams. It is interesting to compare this figure with Fig. 6, where $`^{min}[Q_{a_\nu }(X^i)]`$ is given as a function of time for different metallicities and ions. A similar behaviour is observed in both figures: the larger the difference of the $`^{min}[Q_{a_\nu }(X^i)]`$ values, the lower the correlation coefficient. This is easily understood: a larger difference in $`^{min}`$ means that the stars that contribute the most to the corresponding effective rates of ionising photons are different, and hence, the effective rates are loosely correlated.

In general, all correlation coefficients are close to unity until the WR phase begins. This phase introduces extra ionising photons associated with a handful of luminous stars, and hence destroy the correlations between the effective rates of ionising photons. Correlation coefficients closer to one are preferentially found at low metallicities, due to the paucity of WR stars.

Fig. 8 is divided into five main panels corresponding to different metallicities. Each main panel is further subdivided into four sub-panels: the top-left sub-panel shows $`\rho (\mathrm{H}^0,\mathrm{He}^0)`$, and $`\rho (\mathrm{H}^0,\mathrm{He}^+)`$, the top-right sub-panel shows the correlation coefficients of ions with similar ionising potential, the bottom-left sub-panel shows the correlation coefficients related to the $`\eta `$ parameter, and the bottom-right sub-panel the correlation coefficients related to the \[N ii\] 6584/H$`\alpha `$ vs. \[O iii\] 5007/H$`\beta `$ diagram.

The $`\rho (\mathrm{H}^0,\mathrm{He}^0)`$ and $`\rho (\mathrm{H}^0,\mathrm{He}^+)`$ correlation coefficients have values significatively different from 1 starting from the onset of the WR-phase. The larger the metallicity, the lower the correlation coefficient, implying that the ratios $`Q(\mathrm{He}^0)/Q(\mathrm{H}^0)`$ and $`Q(\mathrm{He}^+)/Q(\mathrm{H}^0)`$ are not robust constraints of the general properties of a cluster (like the upper mass limit in the IMF, the IMF slope, and the age): they are quite dependent on the particular stellar populations in the cluster and their value can strongly vary from cluster to cluster, especially for poorly populated high metallicity regions (see, for example, Bresolin & Kennicutt 2002). On the other hand, for the same reason they can be used to infer the relative abundance of particular populations in a given observed cluster.

The correlation coefficients shown in the top-right sub-panels are $`\rho (\mathrm{H}^0,\mathrm{O}^0)`$ (the two ions having ionising potentials of 13.599 and 13.618 eV respectively), $`\rho (\mathrm{He}^0,\mathrm{S}^+)`$ (ionising potentials of 24.588 and 23.330 eV respectively), $`\rho (\mathrm{He}^+,\mathrm{O}^{++})`$ (ionising potentials of 54.418 and 54.936 eV respectively), and $`\rho (\mathrm{O}^+,\mathrm{S}^{++})`$ (ionising potentials of 35.118 and 34.830 respectively). With the only exception of $`\rho (\mathrm{He}^+,\mathrm{O}^{++})`$, the correlation coefficients are pretty close to 1 in all the cases, as expected during the WR phase. In contrast, the correlation coefficient of $`\rho (\mathrm{He}^+,\mathrm{O}^{++})`$ does not have a constant value, but surprisingly varies between values close to one and values close to zero. The zero value does not necessarily imply that the two quantities are not correlated in the corresponding age ranges: it rather indicates that one or both of the $`Q_{a_\nu }`$ values are zero, that is $`\rho (\mathrm{He}^+,\mathrm{O}^{++})=0`$ implies the absence of photons above 54.418 eV in the model atmospheres of stars at the corresponding ages.

The bottom-left sub-panels correspond to correlation coefficients related to the $`\eta `$ parameter, that is the ratio of $`\mathrm{O}^+/\mathrm{O}^{++}`$ to $`\mathrm{S}^+/\mathrm{S}^{++}`$ (Vílchez & Pagel 1988). Note that the $`\eta `$ parameter is measured from forbidden lines, thus it is related to $`Q_{a_\nu }(\mathrm{O}^0)`$, $`Q_{a_\nu }(\mathrm{O}^+)`$, $`Q_{a_\nu }(\mathrm{S}^0)`$ and $`Q_{a_\nu }(\mathrm{S}^+)`$. In general, the correlation coefficients involving $`\mathrm{S}^0`$ are markedly lower than 1. Remember that $`\mathrm{S}^0`$ was the ion with the lowest $`^{min}`$, which also means that it is produced by a population of stars different from those responsible for the ionising flux (in fact, it is produced by the ionising stars plus other stars): this fact naturally produces a poor correlation. As in the case of $`\rho (\mathrm{H}^0,\mathrm{He}^0)`$ and $`\rho (\mathrm{H}^0,\mathrm{He}^+)`$, the $`\eta `$ parameter is not a good parameter to constrain the age of the cluster, but it is a good parameter to determine the relative proportion of non-ionising and ionising stars, and the softness of the ionising radiation for specific clusters.

Finally, the interesting correlation coefficients $`\rho (\mathrm{H}^0,\mathrm{O}^+)`$ and $`\rho (\mathrm{N}^0,\mathrm{O}^+)`$, together with $`\rho (\mathrm{H}^0,\mathrm{N}^0)`$, can be used to estimate the dispersion in the \[N ii\] 6584/H$`\alpha `$ vs. \[O iii\] 5007/H$`\beta `$ diagram. The effective rates of ionising photons $`Q_{a_\nu }(\mathrm{H}^0)`$ and $`Q_{a_\nu }(\mathrm{N}^0)`$ are strongly correlated, implying that the abundances of the two ions are linearly dependent. In turn, this relation suggests that both elements will show a similar correlation with O+, a fact confirmed by $`\rho (\mathrm{H}^0,\mathrm{O}^+)`$ presenting a behaviour similar to $`\rho (\mathrm{N}^0,\mathrm{O}^+)`$. This has an additional implication: the position of data points in the \[N ii\] 6584/H$`\alpha `$ vs. \[O iii\] 5007/H$`\beta `$ diagram depends on two variables only (either O+ and H0 or O+ and N0), and so the data points will lie on a uniparametric line rather than being spread on the entire plane. Hence, the correlation coefficients obtained through the analysis of the sampling effects may explain the fact that observational data follow a narrow sequence. However, the sequence is uniparametric only for a given metallicity, and it could be different when metallicity varies, yielding a larger scatter in the diagnostic diagram. That is, in fact, the situation in the lower region of this diagram, which is populated by data points corresponding to different metallicities.

## 5 Conclusions

In this paper we have presented an evaluation of the dispersion in the properties of the ionising flux obtained from synthesis models due to the IMF sampling. The ionising fluxes have been characterised by appropriate physical quantities, the effective rates of ionising photons, which are obtained from the integration of the ionising flux weighed by the photoionisation cross section of different ions, with the future goal of using them in the evaluation of sampling effects on emission line spectra.

We have obtained the effective rate of ionising photons for a wide set of metallicities and atmosphere models, comparing the results of two different codes to assess some of the systematic effects. The analysis of the effective rates obtained with different model atmospheres is consistent with previous studies on the impact of the atmosphere models on the ionising flux of star forming regions, and it additionally provides a proper statistical ground for such assessments. We have identified problems in the computation of the ionising flux in regimes where stellar evolution does not present a continuous behaviour with the initial stellar mass, such as at the end of the WR phase.

We have computed the corresponding minimum initial cluster masses for which the ionising continuum obtained by synthesis models can be applied. The minimum masses range from 103 to more than 106 M depending on the metallicity and the age of the stellar population. For observed clusters with stellar masses below the minimum mass, the ionising flux is better accounted for by a single star, rather than by a population of massive stars, so that the use of synthesis models is not appropriate in these cases. We have also evaluated the dispersion in the distribution of the effective rates of ionising photons at the minimum mass, and found it to be around 60%. This means that the initial mass of observed stellar clusters must be larger than a minimum value in the range $`3\times 10^4`$ M to $`3\times 10^7`$ M, depending on age and metallicity, in order for evolutionary synthesis models that assume a completely sampled IMF to be applied. We have also found indications suggesting that the emission of \[S ii\] is much less affected by sampling effects than that of other ions. The He ii lines are in turn the most affected by sampling, especially during the WR phase, and should not be used to constrain the evolutionary status of stellar clusters.

We also studied the correlations between different effective rates of ionising photons, finding in some cases correlation coefficients close to one. This result agrees with the observational finding that H ii regions over a vast range of scales, ranging from those ionised by single stars to those ionised by super-stellar clusters, are found in a relatively narrow band in some emission-line diagnostic diagrams.

Tables containing the effective rate of ionising photons, their corresponding $`𝒩`$, the corresponding LLL, and the complete correlation matrix for H0, He0, He+, N0, N+, N++, O0, O+, O++, S0, S+, S++, C0, Ne0, Ne+, Ar0, Ar+, and Ar++ ions are available at the www address http://laeff.inta.es/users/mcs/SED/. They include the results using the CMHK code with standard and high mass-loss rate444The resulting values in those tables are normalised to 1 M for a Salpeter IMF in the mass range 2–120 M instead the 0.1–120 M mass range used in this paper..

###### Acknowledgements.

We want to acknowledge the referee, L. Smith, for very useful comments. We want also acknowledge G. Ferland for his clear and extensive documentation of Cloudy. MC has been partially supported by the Project AYA3939-C03-01 of the Spanish Programa Nacional de Astronomía y Astrofísica of the MCyT. VL is supported by a Marie Curie Fellowship of the European Community programme Improving Human Research Potential and the Socio-Economic Knowledge Base under contract number HPMF-CT-2000-00949.