{'fraction_non_alphanumeric': 0.02760103181427343, 'fraction_numerical': 0.010346804241903124, 'mean_word_length': 4.309846294323543, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
R-matrix propagation with adiabatic bases for the photoionization spectra of atoms in magnetic fields

The photoionization spectrum of an atom in a magnetic field is calculated by combining R-matrix propagation with local adiabatic basis expansions. This approach considerably increases the speed and the energy range over which calculations can be performed compared to previous methods, allowing one to obtain accurate partial and total cross sections over an extended energy range for an arbitrary magnetic field strength. In addition, the cross sections for all atoms of interest can be calculated simultaneously in a single calculation. Multichannel quantum defect theory allows for a detailed analysis of the resonance structure in the continuum. Calculated cross sections for a range of atoms in both laboratory and astrophysical field strengths are presented.


I. INTRODUCTION
The spectrum of an atom in a magnetic field has played an important role in the development of quantum theory and in atomic structure [1]. For a long time the effect of the external field on the atom was treated perturbatively. However with the discovery of large magnetic fields in white dwarf stars ( 10 2 − 10 5 T ) and neutron stars ( 10 7 − 10 9 T ) in the seventies attention has focussed on non-perturbative treatments of the field-atom interaction [2]. For laboratory strength magnetic fields ( 6T ) the potential due to the applied field only becomes comparable to the intrinsic Coulomb potential for an electron in a high Rydberg state or continuum state of the atom. In the eighties and nineties research was focussed on atoms in laboratory strength fields as these systems have an inherent non separability arising from the competing spherical symmetry of the atom and the cylindrical symmetry of the applied field which leads to the classical system exhibiting chaotic behavior [3]. An atom in a field thus provided an experimentally realizable quantum system whose corresponding classical phase space is chaotic, serving as a prototype for studying classical and quantum chaos. This lead to fruitful developments in the theory of quantum chaos [4].
For bound states of atoms in a magnetic field large scale basis set calculations have proved very successful in finding the energy eigenstates and photoabsorbtion spectra of atoms in moderately strong fields [2]. For super strong fields such as those found in neutron stars, the field can modify the atomic structure of the ground state of the atom and different theoretical techniques have to be used. Currently for these cases the energy eigenstates and photoabsorbtion spectra are only known for some low levels of a few light atoms [5].
The positive energy or continuum spectra of an atom in a field proved more challenging particularly in calculating photoionization cross sections at laboratory strength fields. The three main theoretical methods that have been successful at calculating cross sections at both laboratory and astrophysical strength fields are the complex coordinate method of Delande et al [6], the R matrix method of O'Mahony and Mota-Furtado [7,8], and the diabatic by sector method of Watanabe and Komine [9]. A detailed comparison between theory and experiment has been possible due to the high resolution experiments carried out by Iu et al [10] on lithium in a field of about 6T . Although all three approaches have recreated the experimental spectrum of Iu et al over a narrow energy region they are not particularly suited to calculating the photoionization spectrum over a large energy region. We present here a major improvement on previous R matrix methods applied to this problem by using local adiabatic basis states to propagate the R-matrix from low r to the asymptotic region. This leads to a large saving in both the CPU time and computer memory required to perform the calculation. In addition, for a given value of the field strength, the cross sections for all atoms of interest can be calculated in one step without the need for any additional propagations. We demonstrate that by using a combination of R matrix propagation with local adiabatic basis states and multichannel quantum defect theory (MQDT) [11] one has a method that can be used to calculate photoionization cross sections of any atom over very large energy regions and field strengths. MQDT can also be used to analyze the resonance structure in detail. An efficient approach to calculating such cross sections is of importance in many areas where magnetic fields play a role, for example in calculating stellar opacities for magnetic white dwarfs or for evaluating recombination rates for an atom or ion in a magnetic field at low temperatures (e.g. anti-hydrogen) where one has to calculate the cross section over very large energy regions for a given field strength.
In section II we give the theory used to evaluate the photoionization cross section of an atom in an external magnetic field by combining R matrix propagation with local adiabatic basis states and MQDT. In section III all the required details of the computation are given. The results of the calculations are presented in section IV for a variety of atoms at both laboratory and astrophysical field strengths and a concluding section is given in section V. In the Appendix we give details on how to construct the Hamiltonian matrix for the propagation.

II. THEORY
The Hamiltonian for a hydrogen atom in a magnetic field (taken to be in the z direction) in the symmetric gauge can be written using atomic units ( = m = e = 1) as [1], where the magnetic field B is measured in atomic units by β = B/B 0 with B 0 = 4.70108 × 10 5 T . (We neglect the spin as it only produces a uniform shift in the energy scale). The Hamiltonian has two conserved quantities in addition to the energy, namely the z-component of the angular momentum L z and π z the z-parity. The eigensolutions can thus be studied for fixed values of m the azimuthal quantum number and for π z = ±1. The linear Zeeman term in eqn.
(1) thus only adds a uniform shift to the total energy. When a photon excites an electron to the continuum, for a given magnetic field strength, one can in general identify three regions of interaction with the continuum electron [8]. Typically for low r the spherically symmetric Coulomb potential dominates over the cylindrically symmetric diamagnetic term or quadratic term in eqn. (1), at intermediate values of r the Coulomb and magnetic potentials are of comparable strength (the strong mixing region) and at high values of r or asymptotically in r the cylindrical symmetry of the diamagnetic potential predominates. For a general atom (or molecule) one adds a fourth region, the core, where the excited electron interacts with the multi-electron core before emerging into the Coulomb region described by the Hamiltonian in eqn. (1). Exploiting this natural partition in configuration space forms the basis of the R-matrix approach to solving atomic and molecular problems where solutions are sought in each region and then matched together at the boundaries between the regions to form the solution over all space [12]. A novel aspect of the magnetic field problem is having to deal with the change in symmetry from spherical to cylindrical which involves introducing two-dimensional matching procedures [8]. We describe below how the R-matrix is propagated through the regions described above and how the two dimensional matching procedure is implemented to give the reactance matrix and the photoionization cross section.

A. Propagating the R-matrix
An atom in a magnetic field is assumed to be excited from an initial state, either a ground or low lying excited state, by a polarized photon leading to an electron in the continuum with specific values of m and π z . The electron emerges from the first region, the core region, into the second or Coulomb region with energy ǫ where eqn.(1) can be approximated by the field free Hamiltonian because the diamagnetic terms are negligible in comparison with those from the Coulomb potential of the atomic core. (We shall assume here that the field strengths for non-hydrogenic atoms are not large enough to significantly distort the core (i.e. β < 1). In this case a different treatment would be required for the core although the propagation outside the core could still be implemented). Therefore at some radius r = a in the Coulomb region the radial form of the wavefunction, F ǫ l (r), of the continuum electron can be written in terms of a linear combination of the energy normalized regular s and irregular c Coulomb functions [13] in spherical coordinates or a phase shifted Coulomb function giving the general form of the solution as where the A ǫ l are constants to be determined. The quantum defects µ l represent the effects of the non-hydrogenic field free core [11] and can be calculated ab-initio or taken from experiment. (The more general case of a multi-channel wavefunction with a reactance matrix instead of quantum defects in eqn. (2) is straightforward to include). Y lm (θ, φ) are the spherical harmonics. Knowing the phase shifted Coulomb function and its derivative at r = a, the R matrix [12] can be constructed on the outer boundary of this region as since the coefficients A ǫ l cancel out in the above expression. In the third or strong mixing region the effects due to the Coulomb potential of the core and those due to the magnetic field are of a comparable size. This region is defined by a < r < b where the radius b is taken to be large enough such that the Hamiltonian is separable in cylindrical coordinates. The change in symmetry of the potential, from spherical to cylindrical, is therefore completely contained within this region. We wish to propagate the initial R-matrix in eqn. (3) to obtain the R matrix at the outer boundary of this region r = b [14,15]. To do this we divide the region into N radial sectors with radii a → a 1 , a 1 → a 2 , · · · , a N −1 → b. The size of each sector and the number of sectors N are important parameters in the calculation and we show how these are optimally determined later. Within each of the sectors we construct a local adiabatic basis as follows. For the n th sector we take a radius r n a within the sector a n−1 < r n a < a n , which is usually the mid point of the sector, and we diagonalize the fixed r or adiabatic Hamiltonian H ad in a basis set of spherical harmonics such that where the d lλ are constants giving the eigenvalue equation The adiabatic states form a locally optimized basis set for each sector and the adiabatic potential curves can be produced by plotting the eigenvalues obtained from the diagonalization at successive r n a . A typical set of curves is shown in figure 1. At small r the adiabatic eigenfunctions functions are spherical harmonics and the potential exhibits the centrifugal barrier. At large r these curves have the equal energy spacing of Landau states (combined with the −1/r fall off) and the eigenfunctions are localized in the cylindrical coordinates ρ and φ. The change in symmetry happens predominantly around a region of avoided crossings that can clearly be seen in the diagram. The functions φ λ are therefore a very good basis with which to represent the angular part of the wavefunction in the local region around r n a , namely within the sector n. To propagate the R-matrix from sector to sector we use these basis states within each sector and we diagonalize the full Hamiltonian (eqn. (1)) plus the Bloch operator L or surface term [16], in a basis set consisting of a product of orthogonal radial functions f j (r) and the adiabatic functions generated for that sector φ λ (r n a ; θ, φ). The radial basis functions used are defined in terms of Legendre polynomials P j as follows [15] f j (r) = 2j − 1 a n − a n−1 where u = 2 a n − a n−1 r − a n + a n−1 2 .
The eigenvalue equation is thus The eigenfunctions obtained from diagonalizing this operator are therefore The total continuum wavefunction at any energy ǫ, Ψ ǫ , can be expanded in terms of these R-matrix eigenstates Ψ k . Since the general solution in the n th sector can also be written as Ψ ǫ = λ F λ (r)φ λ (r n a ; θ, φ) it is straightforward to show [15], using the operator H + L and eqn. (9), that the values of the functions and derivatives on the boundaries of the sector are related by F (a n−1 ) = r 2 F ′ (a n ) − r 1 F ′ (a n−1 ) F (a n ) = r 4 F ′ (a n ) − r 3 F ′ (a n−1 ), (11) where the matrix elements of r 1 to r 4 , called the sector R-matrices, are given by and In summary, knowing the eigenvalues and eigenvectors of eqn. (9) one can construct the sector R-matrices r 1 to r 4 above which, through eqns. (11), relate the radial solutions and their derivatives on the boundaries of the sector. The R-matrix relates the function to its derivative, i.e. F (a n ) = R(a n ) F ′ (a n ), and a simple manipulation of eqn. (11) yields the relationship between the R-matrix on the inner and outer boundaries of the sector R(a n ) = r n 4 − r n 3 r n 1 + R(a n−1 ) where R(a n ) and R(a n−1 ) are represented in the same adiabatic basis set as the sector R matrices.
As the adiabatic basis changes from sector to sector we need finally to change the basis representation of the R matrix. The matrix with elements is thus constructed. One uses this transformation to change the basis representation of R giving where T is the transpose of T. Starting with some initial R matrix it can thus be propagated from sector to sector using eqns. (14) and (16). Using the initial input R-matrix given by eqn. (3) the propagation gives the final R matrix on the outer boundary in the asymptotic region at r = b with the final R matrix being represented in the local adiabatic basis of the last sector.
Although it is possible to propagate the R matrix itself at each of the sector radii as described above, it is more practical and efficient to derive global sector R matrices (R 1 , R 2 , R 3 and R 4 ) relating the first and n th sectors [17] which can be built up sequentially using eqn. (11). One initially generalizes eqn. (11) to relate the values of the functions and derivatives on the boundaries of the first and n th sectors The operator relations for these global sector R matrices, including the change in basis, have been derived by Stechel et al [17] and can be obtained by matching the wavefunction and its derivative on the boundaries between each of the n sectors . They are In these equations R n i are the global sector R matrices for all sectors up to n, r n i are the sector R matrices for sector n and T n−1,n is the transformation matrix between the adiabatic basis used in sectors n − 1 and n. In an analogous way to eqn. (14) the global sector R matrices at the end of the propagation through all N sectors are used to relate the R matrix at r = a to the R matrix at r = b Note that R N 1 to R N 4 , determined from eqn. (18) above, are independent of which atom one uses so that once they are calculated from the propagation then R(b) for a whole set of atoms can be evaluated at once using eqn. (19) by just using the appropriate quantum defects to calculate R(a) in eqn. (3). Hence for given values of B, m and π z , the cross sections of all atoms of interest can be calculated simultaneously without the need for any additional propagations.

B. Asymptotic region
Having found the R-matrix at r = b by propagation we must match this to the asymptotic solutions at r = b to find the solution over all space. For large r, r > b, the magnetic field dominates. Since the motion in ρ is bounded − 1 r → − 1 z , and the Hamiltonian in eqn. (1) is separable in cylindrical coordinates where H L is the Hamiltonian for the Landau states H L has eigenvalues E L i = (2i + |m| + m + 1)β, i = 0, 1, 2 . . . and its eigenfunctions are the Landau states Φ i (ρ, φ) [1]. The asymptotic region c ≤ z ≤ ∞, 0 ≤ ρ ≤ ∞ with c less than the radius r = b, is therefore chosen to conform with the cylindrical symmetry of the problem. In this region a set of j linearly independent solutions may be written, as is standard in scattering theory [11], in terms of the solutions of eqn. (20). These are a product of Landau states Φ i and a linear combination of energy normalized regular and irregular Coulomb functions in z, s and c, evaluated at an The constants K kj , the reactance matrix or K matrix, are to be determined by the matching procedure. The R matrix R(b) having being evaluated using eqn. (19) is now matched through a two dimensional matching procedure to these asymptotic solutions on an arc at r = b [8]. To perform this matching the integrals must be evaluated, i.e. the asymptotic solutions must be projected onto the local adiabatic solutions on the radial arc at r = b. This is done by evaluating numerically the four matrices P, Q, P ′ , Q ′ with elements where ′ indicates the derivative with respect to z. This gives the regular and irregular components of the solutions at r = b and allows one to calculate the outer R-matrix from the asymptotic region at r = b in terms of K the reactance matrix. Equating the inner and outer R-matrices at r = b gives Re-arranging this expression gives us the equation for the K matrix Knowing K we have the energy normalized solution over all space and we can calculate both partial and total photoionization cross sections.

Multichannel quantum defect theory
For a given total energy ǫ there are two possibilities for the behaviour of the Coulomb solutions in eqn. (22). If ǫ i = ǫ − E L i > 0 the channel is open and the Coulomb functions oscillate at r = b and all the way to infinity. If ǫ i < 0 then the channel is closed and the solutions must decay as z → ∞. The physical K matrix therefore is a square matrix with the dimension of the number of open channels. However with multichannel quantum defect theory one exploits the known analytic properties of the Coulomb functions to enforce the boundary conditions for the closed channels. For the closed channels there are two possible scenarios at r = b. Either the channel is strongly closed and is already exponentially small at r = b or it is weakly closed, i.e. the Coulomb functions are still oscillating at r = b before decaying at infinity. In MQDT, for a weakly closed channel, one instead uses Coulomb functions in eqn. (22) which don't decay at infinity and one treats the channel as if it is open only enforcing the boundary condition at infinity analytically in a final step [11]. Therefore the resonance structure due to these weakly closed channels can be calculated analytically. Hence the K matrix calculated by doing the matching using MQDT, denoted by K, has dimension given by the number of open plus weakly closed channels. In general K has a smooth dependence on energy because the Rydberg series of resonances converging on the Landau thresholds corresponding to the weakly closed channels have not as yet been included.
The open part of the actual physical reactance matrix K can be recovered from the matrix K by the formula [11] K oo = K oo − K oc tan(πν) + K cc The o and c subscripts refer to the open and closed channels of K and the tan πν ′ s form a diagonal matrix where the ν's are related to the energies ǫ i by This way of obtaining the K matrix has two major advantages. Firstly it is much quicker computationally since K represents a reactance matrix with a lot of the resonance structure removed it varies much more slowly with energy than the full reactance matrix K. This allows one to calculate K on a coarse energy mesh with fairly large energy spacings. These K's can then be used to calculate K over an arbitrarily fine energy mesh using the analytic formula in eqn. (27). This allows the propagation stage to be performed at fewer energy points ultimately speeding up the calculation enormously. The second major advantage of this approach is that the resonance structure converging to a particular Landau threshold can be identified by removing it from K. This is done by keeping open the relevant Landau channel in the evaluation of equation (27). By comparing a spectrum with all resonances converging on a particular Landau threshold removed with that of the full spectrum it is possible to determine which resonances converge to which thresholds [8,18]. This technique is demonstrated in section IV.

C. Photoionization cross section
The photoionization cross section is given by where α is the fine structure constant, ω the photon energy, ǫ the polarization direction, Ψ o the initial bound state and Ψ − ǫ the energy normalized 'incoming' wavefunction [19]. There is a standard transformation to go from the K-matrix form given in (22) to the S-matrix or incoming form Ψ − ǫ [11]. Once the K matrix and hence the asymptotic form is known from the matching, the photoionization cross section is evaluated by calculating the amplitude of the wavefunction near the origin, and hence the dipole integrals, as follows.
The 'incoming' wavefunction on the inner boundary of the strong mixing region (r = a) can be written as The wavefunction on the outer boundary of the strong mixing region can be written as where G − λ is the energy normalized 'incoming' asymptotic solution calculated from the matching procedure in section II B. The radial solution F − and G − are linked by the global sector R matrices as given in eqn. (17) The radial solutions on r = a can be written as F − = S A − from eqn. (2) where S is a diagonal matrix with elements s ǫ l (a) + c ǫ l (a) tan πµ l and A − are the field and energy dependent amplitudes. Substitution into equation eqn. (32) yields the equation The coefficients A − can therefore be evaluated by solving the set of linear equations because G ′− (b) is known once the K matrix is known. Using the coefficients A − one can then simply express the cross section in terms of the field free photoionization cross section as the coefficients A − give the difference in amplitudes between the field dependent and field free amplitudes . For example for excitation from the 1s state of hydrogen using linear polarized light (∆m = 0) only the l = 1 component of the final state will be accessed. The photoionization cross section in this case is therefore given by where σ B=0 is the field free photoionization cross section for hydrogen.

III. COMPUTATIONAL DETAILS
The first point to address is the choice of the radii a and b. The inner radius must be taken larger than the core (∼ 1 a.u.) for non-hydrogenic atoms yet small enough that the diamagnetic term is still negligible compared to −1/r. For example we took a = 200 atomic units for laboratory strength fields. For the asymptotic radius b, extensive calculations show that it is surprisingly large indicating that the long range coupling due to the terms of O(1/z 3 ) falls off very slowly. Based on our experience we used the following empirical formula for an arbitrary field b = 700 The number of sectors between r = a and r = b and their sizes are chosen in the following way. First one chooses the maximum energy to be used in the calculation, ǫ max . The transformation matrix (T n−1,n ) λλ ′ between adiabatic functions evaluated at two radii r n−1 a and r n a for consecutive sectors can be calculated using eqn. (15). If r n−1 a was equal to r n a , (T n−1,n ) λλ ′ would be the identity matrix. As the distance | r n−1 a − r n a | increases the corresponding adiabatic functions become more different. The effect this has on T is that the diagonal elements get smaller and the off diagonal elements get bigger. The maximum sector size for the angular adiabatic basis is found by choosing a limit on how small any of the diagonal elements of T can become. Choosing a value of 0.5 for the smallest diagonal element restricts the size of each sector for a laboratory strength field to those shown in figure 2. One can see that for small and large r the adiabatic functions don't vary much with r being close to spherical harmonics and Landau states respectively leading to large sector sizes. The intermediate range of r, where the avoided crossings in the potential curves shown in figure 1 are present, is where the angular functions are changing rapidly and one requires small sectors.
The radial basis sector widths are chosen by comparison with the local wavelength of the Coulomb functions, namely Radial sector size = constant If a radial basis of ten Legendre polynomials is taken, it is found that a constant = 6.0 is sufficient to produce an accurate description of the wavefunction [15]. For the maximum energy used, the sectors resulting from this criterion are also shown in figure 2. There are therefore two criterion for choosing the sector sizes; one from the radial basis and one from the angular basis. The smallest value of these two determines the actual sector sizes and these are also shown in figure 2. The adiabatic functions φ λ (r n a ; θ, φ) are obtained by diagonalizing the adiabatic Hamiltonian H ad in a basis of spherical harmonics. The number of spherical harmonics used must ensure the functions φ λ (r n a ; θ, φ) are properly converged. At the matching radius r = b, the ratio between the diamagnetic term and the potential term in H ad will be at its greatest. It is therefore at this radius that the largest number of spherical harmonics will be required because it is here that the adiabatic functions will have their greatest degree of cylindrical symmetry. We choose the number of l's so that the eigenvalue corresponding to the second closed channel at r = b for an energy ǫ max is accurate to 0.5%. It has been verified that this is sufficient to give convergence in the final cross section.
When the R matrix is propagated from r = a → r = b, it is necessary to retain, within any one sector, all of the locally open channels plus a few of the locally closed ones. The number of closed channels required is directly related to the threshold value on the adiabatic angular functions used to determine the sector sizes. The number of closed channels needed, however, is constant for all of the sectors. For a given energy ǫ, because the number of locally open channels changes with radius (see figure 1), the total number of channels retained in any one sector varies. Because it is very difficult to include a channel half way through the propagation, the number of channels retained in a sector is found by the following method.   The number of channels needed to match the R matrix to asymptotic solutions in order to obtain the K matrix was checked for each spectrum to ensure convergence. Once the K matrix was obtained, however, the number of channels needed for the calculation of the cross section was reduced to the number of open channels plus the number of weakly closed channels. Weakly closed channels were taken to be those that were within an energy (in a.u.) of 1 b of their corresponding Landau energy.
The cross section is obtained over a coarse mesh of energies initially and then MQDT is used to calculate the cross section over an arbitrarily small energy range. The convergence of the cross sections was tested by varying the inner and outer radii , a and b, and the other parameters in the calculations.

IV. PHOTOIONIZATION CROSS SECTIONS FOR LABORATORY STRENGTH FIELDS
Calculations at laboratory strength fields are the most demanding numerically as the radial distances over which one needs to propagate the R-matrix are large. The number of channels needed in the calculation can also become quite large. We focus on the photionization spectrum of lithium in a magnetic field of 6.1143T (β = 1.3 × 10 −5 ) where experimental data exist (Iu et al [10]). Lithium is excited from the 3s state using linearly polarized light giving a final state with m = 0 and π z = −1 or odd z-parity. We show in figure 3 the calculated spectrum between the photoionization threshold (i.e. i = 0 or first Landau level) and the second Landau level or threshold i = 1.
The only significant quantum defect is for l = 1 , µ p = 0.053. (The l = 0 continuum state doesn't play any role for excitation from s states.) The photoionization cross section is given in arbitrary units but it can be put on an absolute scale if the field free cross section is known. The radius a is 200 and b = 12600. The adiabatic matrix threshold was taken to be 0.1. Together with an ǫ max of 3.9 × 10 −5 a.u. this gives the number of sectors to be 64. The maximum number of locally open adiabatic channels is 27 and taking 13 extra closed channels the maximum number of channels overall is 40, hence the maximum size matrix to be diagonalized is 400 since there are 10 radial basis functions per sector. The cross section is calculated over a course mesh at 400 energy points between the thresholds. The full spectrum is then calculated semi-analytically using MQDT and eqn. (27) with over 10,000 energy points. The spectrum obtained from this calculation is displayed in the top panel of figure 3 showing the complete resonance structure between the thresholds. MQDT can also be used to examine the resonance structure by keeping individual weakly closed channels open, i.e. by not enforcing the closed channel conditions asymptotically. For a single threshold this is equivalent to a Gailitis average [11,20] over the resonances converging to that threshold. For several thresholds this is a type of generalized Gailitis average over the resonance structure and would be similar to convoluting the actual spectrum with a Gaussian of a certain width [18]. The differences between the first and second panel allows In addition to helping in the analysis of the spectrum, MQDT allows one to calculate the resonance structure to an arbitrarily small resolution. This is demonstrated in figure 4 where an enlargement of the spectrum in a very small energy region just below the first excited Landau level is shown.
However the real power of the combination of R matrix propagation with adiabatic bases can be seen in figure 5, where we calculate the photoionization cross section over an extended energy region covering over 6 Landau thresholds from the ionization threshold. Because the number of open channels increases with energy the size of the matrices to be diagonalized are slightly larger. The total number of sectors increases also but even for the highest energy, the number of channels didn't exceed 50 giving the largest matrices that need to be diagonalized to be of the order of 500. The cross section was calculated on a coarse energy grid of 500 points over each of the thresholds before applying MQDT. The method scales in a reasonable way thus allowing one to calculate the photoionization cross section of an atom over a large energy range above the ionization threshold.

V. PHOTOIONIZATION CROSS SECTIONS FOR ASTROPHYSICAL STRENGTH FIELDS
Photoexcitation and photoionization cross sections of light elements such as hydrogen and helium are important in understanding the properties of white dwarf and neutron stars [2]. The method detailed in section II can equally well be applied to atoms in astrophysical strength magnetic fields. In fact the computation is much easier in this case compared to laboratory strength fields as the radius b is much smaller and the number of channels and sectors is smaller too. We give just two examples to illustrate the point. We first consider the photoionization spectrum of hydrogen from the ground state using linearly polarized light in a field of 23,500T (β = 0.05 a.u.) (See figure 6.) The radius a was 1 a.u. and b = 50, the number of sectors used was 20 and up to ten channels were used in the propagation. This spectrum has also been calculated using the complex co-ordinate method [6,21]. Excellent agreement is found  between the two methods away from the energy region near to the ionization thresholds. As the complex co-ordinate method uses a finite basis it cannot represent all of the Rydberg structure just below the ionization thresholds.
The second example is the photoionization spectrum of helium from the ground state in a field of 4700T (β = 0.0001 a.u.). The only significant quantum defect is for l = 1, µ p = −0.012. (The l = 0 continuum state will only play a role for excitation from 1s2p state.) The radius a was 5 a.u. and b = 150, the number of sectors used was 20 and a maximum of twelve channels were used in the propagation. Resonances converging to the first 4 excited Landau thresholds are shown in figure 7. The resonance structure is that of strongly perturbed Rydberg resonances converging to the individual thresholds.
The spectra in both cases, when calculated over an extended energy range, repeat the patterns shown in figures 6 and 7.

VI. CONCLUSIONS
We have presented a detailed description of a method to evaluate the photoionization cross section of an atom in an external magnetic field. By combining R-matrix propagation with local adiabatic bases we have shown that is possible to calculate the cross section over an extended energy range for a range of field strengths. We have calculated cross sections for a range of atoms in both laboratory and astrophysical field strengths to illustrate the generality of the method. In addition, for given values of B, m and π z , the spectra of all atoms of interest can be calculated without the need for any additional propagations. By using MQDT one is able to calculate quickly all of the resonances in the spectrum and to analyze some of their main characteristics. Partial cross sections to individual Landau levels are evaluated when calculating the total cross section enabling one to calculate their distributions. The method can be used to calculate the large amounts of data needed for such problems as stellar opacities or for calculating re-combination rates at low temperatures for an atom in a magnetic field.

VII. APPENDIX
In constructing the matrix representation of the Hamiltonian plus Bloch operator in each sector one has to construct matrices for the operators and r 2 with the radial basis set of shifted Legendre polynomials.
For an arbitrary sector with inner radius r = a and outer radius r = b the coordinate r can be rescaled to the coordinate u such that The limits of the integrals in r over the limited range a → b therefore become −1 → 1 in the coordinate u. The orthogonal radial functions are thus where P n is the Legendre polynomial of order n. The kinetic term can be evaluated by integration by parts and using the relation (Copson [22] page 282) where m ≥ n.
The first values could also be calculated using hyper-geometric functions via Q n (c) = π for | c |> 1, however in the calculations described in this paper the expression in eqn. (56) was used. The Legendre functions of the second kind can therefore be evaluated using the recurrence relation Q n (c) = n + 2 n + 1 Q n+2 (c) − 2n + 3 n + 1 cQ n+1 (c). (59)