J Kovář

O Kopáček

V Karas

Z Stuchlík

Off-equatorial orbits in strong gravitational fields near compact objects -II. Halo motion around magnetic compact stars and magnetized black holes

Keywords: numbers: 04.25.-g, 04.40.Nr, 04.70.Bw, 04.90.+e

Off-equatorial circular orbits with constant latitudes (halo orbits) of electrically charged particles exist near compact objects. In the previous paper, we discussed this kind of motion and demonstrated the existence of minima of the twodimensional effective potential which correspond to the stable halo orbits.

Here, we relax previous assumptions of the pseudo-Newtonian approach for the gravitational field of the central body and study properties of the halo orbits in detail. Within the general relativistic approach, we carry out our calculations in two cases. Firstly, we examine the case of a rotating magnetic compact star. Assuming that the magnetic field axis and the rotation axis are aligned with each other, we study the orientation of motion along the stable halo orbits. In the poloidal plane, we also discuss shapes of the related effective potential halo lobes where the general offequatorial motion can be bound. Then we focus on the halo orbits near a Kerr black hole immersed in an asymptotically uniform magnetic field of external origin.

We demonstrate that, in both the cases considered, the lobes exhibit two different regimes, namely, one where completely disjoint lobes occur symmetrically above and below the equatorial plane, and another where the lobes are joined across the plane. A possible application of the model concerns the structure of putative circumpulsar discs consisting of dust particles. We suggest that the particles can acquire a small (but non-zero) net electric charge, and this drives them to form the halo lobes.

Introduction

Motion of matter near compact stars and black holes has been discussed thoroughly in textbooks [START_REF] Chandrasekhar | The Mathematical Theory of Black Holes[END_REF][START_REF] Misner | Gravitation[END_REF]. In the test-particle approximation, the hydrodynamical terms are neglected, as expected in circumstances when the medium is rarefied and the mean free path is comparable with the typical length-scale of the system [START_REF] Prasanna | [END_REF][4][5]. The interplay between gravitational attraction and the electromagnetic attraction or repulsion is essential for characteristics of the motion, namely, its stability properties. Motivation for these studies arises from the problem of motion and acceleration of matter (charged particles or dust grains) [6][7][8][9][10].

Compact stars can be imagined as being endowed with magnetic dipoles that are anchored in the stars and co-rotate with them [6,[START_REF] Mestel | Stellar Magnetism[END_REF][START_REF] Rezzolla | [END_REF]. Rotation gives rise to the electric component of the field, and we assume that the resulting structure dominates over the small-scale turbulent fields. On the other hand, uncharged black holes do not support their own intrinsic magnetic fields; they can only be embedded in external fields [13,14]. In general, the gravitational and electromagnetic fields interact with each other and the result is determined through the Einstein-Maxwell equations. To avoid such a complex description, one can also consider a test electromagnetic field influenced by the given gravitational field.

In both cases, magnetic compact stars as well as magnetized black holes ‡, the interplay of strong gravity, rotation and magnetic field is capable of accelerating electrically charged particles and, possibly, it can also introduce new families of stable motion. The latter will be the main subject of the present paper.

Recently, we have studied the existence of stable circular orbits with fixed latitudes of charged particles moving off the equatorial plane (halo orbits) [START_REF] Kovář | Classical and Quantum[END_REF][START_REF] Stuchlík | Proc. of the Int. Astr. Union, IAU Symposium[END_REF]. We assumed an axially symmetric (aligned) magnetic dipolar field and built our discussion on classical studies [START_REF] Howard | [END_REF]18] derived in the context of weak (Newtonian) gravitational fields, and designed for planetary studies. Particles on halo orbits are bound to the central object by a combined effect of gravitational and electromagnetic forces and they do not have to cross neither the equatorial plane nor the rotation axis. In the poloidal plane, the halo orbits are located in a kind of lobes that are positioned symmetrically above and below the equatorial plane. The lobes may be disjoint or they may be interconnected, depending on the parameters. Our study [START_REF] Kovář | Classical and Quantum[END_REF], within the pseudo-Newtonian approach [19][20][21], has shown that the stable halo orbits can indeed emerge also near magnetic compact stars, where their structure is affected by strong gravitational field. We have also examined the case of charged and rotating (Kerr-Newman) black holes and naked singularities. Surprisingly, the stable halo orbits do not exist above the horizon of Kerr-Newman black holes. This suggests that the structure of the Kerr-Newman spacetime is in several aspects quite special one, while the generic case of dipole-type field usually allows for the motion along stable halo orbits. ‡ We use the term magnetic compact stars (magnetic stars) to specify neutron stars, Q-stars and hybrid stars as well, and the magnetized black holes for the black holes embedded in ordered magnetic fields.

Here, we extend our investigation of the halo orbits to the case of magnetic compact stars and magnetized black holes in general relativity. To this end, we study the existence and orientation of stable halo orbits which surround the minima of two-dimensional effective potential in the poloidal (r, θ)-plane. We also deal with qualitatively different cases of the general, in the halo lobes bounded, off-equatorial motion (halo motion). Namely, we point out to the motion in two separated halo lobes extending symmetrically above and below the equatorial plane versus the motion in a single region arising from the merged lobes across the equatorial plane. We also distinguish the cases when the halo lobes admit outflow of particles to the outer space or they allow the inflow onto the central body.

In the space environment, dust particles can carry small electrostatic charges [22]. An astrophysical motivation to study the halo motion of such weakly charged particles is the search for circumpulsar debris that could be formed from a fall-back disc resulting from the fraction of the explosion ejecta material that fails to escape [23][24][25]. The anomalous X-ray pulsars [26] represent a category of magnetars around which the infrared observation indicate the presence of debris dust discs, potentially relevant for our deliberation. These are young neutron stars in the category of magnetars with slow spin periods in a range ≈ 2-12 seconds. The minimum distance to which the dust can reach is uncertain, however, part of it may enter the magnetosphere and influence the current flows below the light cylinder, until the dust evaporates. The potential relevance of the halo motion of dust particles is based on the fact that these orbits occupy areas of the stable motion. The irradiated dust particles acquire a small net electric charge and modulate the source radiation at their characteristic oscillation frequencies as they move in the halo lobes.

The paper is organized as follows. In section 2, we combine the standard formulation based on the construction of the super-Hamiltonian and the corresponding effective potential, with perhaps a less frequent (in this context) formalism of forces in the projected three-space [27][28][29][30]. We maintain the assumption about axial symmetry but we generalise ref. [START_REF] Kovář | Classical and Quantum[END_REF] by considering the background of the magnetic compact star in a more consistent framework, i.e. in the Schwarzschild geometry with a test dipoletype rotating magnetic field [6] (section 3). Furthermore, as a qualitatively different situation, we also consider a rotating black hole immersed in the uniform magnetic field [14] (section 4). We discuss the astrophysical relevance in section 5 and then we conclude the paper in section 6.

Formalism

According to the standard approach to a test particle motion, we start by construction of the super-Hamiltonian [START_REF] Misner | Gravitation[END_REF]

H = 1 2 g ij π i -qA i π j -qA j , (1)
where m and q are the rest mass and electric charge of the particle, π i is the canonical momentum, and A i denotes the vector potential related to the electromagnetic tensor by F ij = A j,i -A i,j . The particle motion is governed by Hamilton's equations

dx i dλ = ∂H ∂π i , dπ i dλ = - ∂H ∂x i , (2)
where λ = τ /m is the affine parameter and τ is the proper time. § The first Hamilton's equation of motion implies

p i ≡ dx i dλ = π i -qA i . (3)
The second Hamilton's equation ensures that the generalized momenta,

π t = p t + qA t ≡ -Ẽ, π φ = p φ + qA φ ≡ L, (4)
are constants of motion, reflecting the stationarity and axial symmetry of the system. These are connected with Killing vector fields η i = δ i t and ξ i = δ i φ . We start by writing the normalization condition, m 2 = -g ij p i p j , and defining the specific energy, angular momentum and charge, E = Ẽ/m, L = L/m and q = q/m, respectively. We find the two-dimensional effective potential for the particle motion in the form [START_REF] Misner | Gravitation[END_REF]

V eff = -β + (β 2 -4αγ) 1/2 2α , (5)
where

α = -g tt , β = 2 g tφ L -qA φ -g tt qA t , (6)
γ = -g φφ L -qA φ 2 -g tt q 2 A 2 t + 2g tφ qA t L -qA φ -1, (7)
reflecting the motion properties.

Alternatively to the Hamilton's equations [START_REF] Misner | Gravitation[END_REF], one can describe the motion by the Lorentz equation

u k ∇ k u i = q F i k u k . (8)
For our purposes, we find this equation particularly well suited, when being rewritten in the formalism of forces [27]. The forces formalism is based on the projection h ik = g ik + n i n k of the Lorentz equation ( 8) onto the three-dimensional hypersurface orthogonal to the four-velocity field of the Locally Non-Rotating Frames (LNRF) [31],

n i = e -Φ (η i + Ω LNRF ξ i ), e 2Φ = -(η i + Ω LNRF ξ i )(η i + Ω LNRF ξ i ), (9)
where the angular velocity Ω LNRF = -g tφ /g φφ . In the case of a static spacetime, Ω LNRF = 0 and the LNRF become static ones. § We use the geometric system of units (c = G = 1) and a positive signature of the metric. In order to reduce the number of parameters in our classification, we scale all quantities that have the dimension of the power of length by the mass M * of the central object. In this way we adopt the scaled dimensionless quantities. Thus, our formulae become completely dimensionless. Furthermore, our attention is paid to the stationary and axially symmetric spacetimes, described by using the standard Boyer-Lindquist coordinates x i = (t, φ, r, θ), and endowed with the magnetic fields which adopt the same symmetries.

Let us consider an example of the halo motion -the circular motion at constant latitude θ (outside the equatorial plane). The four-velocity field of particles uniformly circling along the halo orbits can be decomposed as

u i h = γ(n i + v h τ i ). (10)
Here, γ = (1 -v 2 h ) -1/2 is the Lorentz factor, τ i = ξ i g -1/2 φφ is the unit spacelike vector orthogonal to n i , along which the spatial velocity v i h = v h τ i is aligned. Both vectors in the decomposition (10) correspond to the base vectors of the standard orthonormal tetrad attached to LNRF: n i = e i (t) and τ i = e i (φ) . Thus, v h is the orbital (azimuthal) velocity measured with respect to LNRF.

Projection of the Lorentz equation, h k j u i ∇ i u k = qh i j F ik u k , can be written in the form

G j + (γv h ) 2 Z j + γ 2 v h C j = -qγ(E j + v h M j ), (11)
where the so-called mass and velocity independent parts of the gravitational, centrifugal and Coriolis inertial forces, and the charge and velocity independent parts of the electric and magnetic forces can be expressed as

G j = -∂ j Φ, (12)
Z j = 1 2 g -1 φφ e -2Φ e 2Φ ∂ j g φφ -g φφ ∂ j e 2Φ , (13)
C j = g -3/2 φφ e -Φ g φφ ∂ j g tφ -g tφ ∂ j g φφ , (14)
E j = e -Φ Ω LNRF ∂ j A φ + ∂ j A t , (15)
M j = g -1/2 φφ ∂ j A φ , (16)
where only the radial and latitudinal components are nonzero.

Our investigation of the existence and orientation of the stable halo motion is based on the following approach. First, we notice that loci of the halo orbits (stable as well as unstable) correspond to the stationary points of the effective potential (5). We are interested in the stable motion mainly, i.e. in those orbits which satisfy the conditions for the local minima of the potential,

∂ 2 r V eff (r, θ; p, L h , q h ) > 0, ( 17
) det H(r, θ; p, L h , q h ) > 0. ( 18
)
Here, H is the Hessian matrix, p is a parameter characterizing the compact object (angular velocity Ω in the case of magnetic star, spin parameter a in the case of Kerr black hole), and L h and q h are the specific angular momentum and charge which are characteristic for a certain halo orbit.

According to the forces formalism, we can determine L h and q h from the balance equations

G r + (γv h ) 2 Z r + γ 2 v h C r = -q h γ(E r + v h M r ), (19)
Note that alternatively to equations ( 19)- (20), the specific angular momentum L h and charge q h can be determined also from the conditions ∂ r V eff = 0 and ∂ θ V eff = 0. However, this routine is more complicated, especially in our second case of the Kerr black hole in the uniform magnetic field.

G θ + (γv h ) 2 Z θ + γ 2 v h C θ = -q h γ(E θ + v h M θ ). ( 20
)
Eliminating q h and assuming 0 < θ < π, θ = π/2, we get the cubic equation

Av 3 h + Bv 2 h + Cv h + D = 0, (21)
where

A = M θ (G r -Z r ) + M r (Z θ -G θ ), (22)
B = E r (G r -Z r ) + E r (Z θ -G θ ) + C θ M r -C r M θ , (23)
C = C θ E r -C r E θ + G θ M r -G r M θ , (24)
D = E r G θ -E θ G r . (25)
The cubic equation (21) has in general three complex solutions,

v h,i = v h,i (r, θ; p), (26)
where i ∈ {I, II, III}. These can represent values of orbital velocities of the charged particles moving along the halo orbits. The corresponding specific charges can be derived from one of the equations ( 19)- (20). We find

q h,i = G r (v 2 h,i -1) -v h,i (C r + v h,i Z r ) (E r + v h,i M r ) 1 -v 2 h,i 1/2 . ( 27
)
Finally, the three values of the specific angular momentum are

L h,i = γ v h,i g 1/2 φφ + q h,i A φ . (28)
Searching systematically through the parameter space (r × θ × p), we can now investigate the validity of conditions [START_REF] Howard | [END_REF], (18) and v h,i ∈ R, |v h,i | < 1, necessary and sufficient for the existence of the stable halo orbits.

Magnetic star

To a good approximation, the gravitational field outside a compact star is described by Schwarzschild metric [START_REF] Misner | Gravitation[END_REF]

ds 2 = -1 -2r -1 dt 2 + 1 -2r -1 -1 dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ). (29)
Naturally, the gravitational field of rotating compact star differs from the Schwarzschild metric (in the case of slow rotation it is well determined by the Hartle-Thorne metric [32]). However, one does not need to consider this discrepancy, taking advantage of the simple analytical form of the Schwarzschild metric element that captures the essential properties of the motion not only near non-rotating black holes, but above the surface of ultra-compact stars and farther away from slowly rotating black holes as well.

The test dipole magnetic field rotating with angular velocity Ω in the Schwarzschild geometry can be expressed in terms of the vector potential [6]

A t = 3 8 ΩMR sin 2 θ, A φ = -3 8 MR sin 2 θ, (30)
where

R = 2 + 2r + r 2 ln 1 -2r -1 . (31)
The related dipole magnetic moment is [33]

M = 4R 3/2 (R -2) 1/2 6(R -1) + 3R (R -2) ln (1 -2R -1 ) B 0 . (32)
This configuration can be considered as a model of magnetic compact star (radius R and angular velocity Ω) endowed with an aligned co-rotating magnetic field of strength B 0 . The latter is measured with respect to the standard orthonormal tetrad of static observers on the surface of the star in the equatorial plane. Rotation terms are neglected in the metric and the frozen-in condition for the magnetic field is imposed, F i j u j MF = 0 (force-free approximation), where 29) and ( 30), the effective potential (5) adopts the form

u i MF = (u t MF , u φ MF , 0, 0), i.e. u φ MF /u t MF = Ω. Employing equations (
V eff = -3 8 qMΩR sin 2 θ + 1 -2r -1 1/2 1 + L r sin θ + 3qMR sin θ 8r 2 1/2 . ( 33
)
The only restriction following from the potential formula is r > 2 (event horizon). However, the effective potential is not meaningful outside the region between the considered surface of the star and the light cylinder, arising due to the condition g ik u i MF u k MF < 0, thus being implicitly defined by the formula

Ω 2 r 3 lc sin 2 θ lc -r lc + 2 = 0. ( 34
)
Table 1. Signs of quantities characterizing the motion along stable halo orbits. Values of r and Ω are taken from the specific regions corresponding to figure 1.

Reg. v h,I Q h,I L h,I A + - + B + - + Reg. v h,II Q h,II L h,II A - + - B - + - Reg. v h,III Q h,III L h,III C + + + C - - -

Existence and orientation of halo orbits

The general analysis in section 2 demonstrated that the dipole magnetic field (30) in the Schwarzschild geometry (29) allows for the existence of the stable halo orbits. We took advantage of the fact that q appears in all formulae as a product with M, so we could introduce the effective specific charge Q = qM instead of the usual specific charge q. This decreases the number of parameters. Only the parameter p ≡ Ω remains in the relations for v h,i (26), L h,i (28) and Q h,i , being expressed by the same formula as q h,i

.

The following conclusions about the existence of stable halo orbits arise from figure 1. The stable halo orbits are possible for particles which co-rotate with the dipole and whose charge satisfies either the condition Q h Ω < 0 (regions A, B ) or Q h Ω > 0 (regions C, C ). By co-rotation of particles we mean v h Ω > 0, i.e. sgn(v h ) = sgn(Ω). Particles can also counter-rotate with the dipole along stable halo orbits; in such case The range and orientation of particle motion along stable halo orbits near a magnetic star. The star rotates with angular velocity Ω. This figure has been constructed for the latitude θ = π/3 and it represents a typical result for the effective potential V eff (r, θ; L, Q, Ω), based on the behaviour of three orbital velocities v h,I (r, θ; Ω), v h,II (r, θ; Ω) and v h,III (r, θ; Ω) (roots of equation ( 21)). These correspond to the specific angular momenta L h,I (r, θ; Ω), L h,II (r, θ; Ω), L h,III (r, θ; Ω) and the effective specific charges Q h,I (r, θ; Ω), Q h,II (r, θ; Ω), Q h,III (r, θ; Ω). The gray areas correspond to the parameters for which the effective potentials V eff (r, θ; L h,I , Q h,I , Ω) (left panel), V eff (r, θ; L h,II , Q h,II , Ω) (middle panel) and V eff (r, θ; L h,III , Q h,III , Ω) (right panel) develop local minima. Only the regions under the light cylinder are considered. The thick vertical lines denote position of the event horizon at r = 2. Description of the regions is summarized in table 1.

the charge must satisfy the condition Q h Ω > 0 (regions A , B). Thus, in the case of the same orientation of the magnetic dipole and its rotation (MΩ > 0), negatively charged particles moving along stable halo orbits can only co-rotate, whereas positively charged particles can either co-rotate or counter-rotate with the dipole. We point out that the above-mentioned results are general. Qualitatively the same structure appears at any latitude, 0 < θ < π, θ = π/2 (the assumption of θ = π/3 in the figure 1 is not crucial). Note that v h L h > 0 holds always for every stable halo orbit.

Effective potential and classification of halo motion

As we can conclude from equation (33), the effective potential is invariant under combinations of simultaneous sign reversals maintaining

sgn(Q Ω) = const, sgn(ΩL) = const. ( 35
)
At the event horizon, the potential diverges as

lim r→+2 V eff = sgn(Q Ω)∞. ( 36
)
In the non-rotating case, the effective potential is invariant under the signs reversals which satisfy

sgn(QL) = const. ( 37
)
The effective potential (33) exhibits a rich spectrum of behaviour. Different topology of the effective potential determines different possible regimes of charged particles motion and yields the following classification. Note that locations of the surface of the star and of the light cylinder must be incorporated into the classification as well. For instance, with respect to the different positions of the light cylinder, we obtain an additional type IIIa to the type IIIa (although the behaviour of the potential itself is identical in both the cases). Similarly, different positions of the surface of the star bring the additional type Ia to the type Ia, and the type IIb to the type IIb. On the other hand, the different types of the potential behaviour itself do not have to enrich the considered classification necessarily. For example, the different number of saddle points that occur below the star surface or beyond the light cylinder does not influence the topology between surface of the star and light cylinder, and is irrelevant for us.

The family of different types of the effective potential can be parameterized by (conserved) energy of the moving particle. Starting at the energy corresponding to local minima of the potential and rising gradually the energy level (and keeping the other parameters constant), we first observe the gradually growing lobes of halo orbits, where particles are trapped. The following classification emerges (see figures 2-3). ¶ Type Ia. Halo lobes merge together once the energy level of the saddle point in the equatorial plane is reached; particles moving within these merged lobes can cross the equatorial plane. Later, when the energy is high enough, particles start falling onto the surface of the magnetic star in the regions where the merged potential lobes open, and the particles hit the surface of the star. Increasing the energy further, the merged lobes open out also on the outer side, crossing the light cylinder. This way the particles can scatter away from the system. Type Ib. Halo lobes merge through saddle points out of the equatorial plane. When the energy is high enough, the merged lobes open out through the surface of the star, and then also through the light cylinder.

Type IIa. Halo lobes merge through the saddle point in the equatorial plane. Then again, the merged halo lobes open out through the light cylinder and, finally, through the surface of the star. Type IIIc. Halo lobes open out through saddle points out of the equatorial plane and merge this way with each other. Later, another corridor between the lobes emerges through the saddle point in the equatorial plane. The potential barrier between the saddle points can be overcome after increasing the energy further. Finally, the merged Type IIb . Halo lobes open out through the surface of the star, and later, when the energy is high enough, through the light cylinder. Increasing the energy further more, the opened lobes merge through the saddle point in the equatorial plane hidden beyond the light cylinder.

Features of the classification

In the case of parallel orientation of the rotation and the magnetic moment of the dipole, the behaviour of the potential of the types Ia, Ib and Ia (serie I) determine the motion of negatively charged particles which co-rotate along the stable halo orbits. At the event horizon, the potential diverges as V eff → -∞, and there are two additional saddle points of the potential near the event horizon (hidden under the surface of the star). In the equatorial plane, the effective potential develops local maxima (also hidden under the surface of the star) in addition to the saddle point. Note that the extent of halo lobes = -5.72 and the corresponding poloidal projection of the motion (lower panels), moving in the separated and merged halo lobes in the vicinity of a magnetic star with R = 4 and Ω . = 0.0115. The trajectories were calculated numerically by using the Lorentz equation of motion (8) with the initial conditions L . = 0.876, φ(0) = 0, r(0) = 5, θ(0) = π/3, and E = 0.8481, u θ (0) = 7 × 10 -4 , u r (0) . = 0.021 (motion in separated halo lobe), and E = 0.8485, u θ (0) = 7 × 10 -3 , u r (0) . = 0.020 (motion in merged lobes).

can never be large (compared to the radius of the black hole horizon) and they cannot be placed farther off the star surface.

Behaviour of the types IIa, IIb and IIb (serie II) determines the motion of positively charged particles, counter-rotating along halo orbits. At the event horizon, the potential diverges as V eff → ∞, and there are no other stationary points hidden under the surface of the star, neither out of the equatorial plane nor within the plane. A relatively high potential barrier develops between the halo lobes and the star. Now, the halo lobes can be of large extent and they can be placed farther away from the star surface.

Finally, the behaviour of types IIIa, IIIb, IIIc and IIIa (serie III) determines the motion of positively charged particles, co-rotating along halo orbits. At the event horizon, the potential diverges as V eff → ∞, and there are no other stationary points out of the equatorial plane. However, in the equatorial plane, the behaviour of the potential is more complicated. In some cases, there can be even more than four stationary points.

There is an infinite potential barrier at the event horizon in the types of series II and III, and behaviour of the potential closer to the magnetic star can be complicated, especially in types of the serie III. But having the sufficiently high energy, particles can always fall onto the star surface, for the star surface always takes place above the event horizon.

The types Ia, IIa, Ia , IIIa and IIIa are characteristic for the 'lower' stable halo orbits, taking place near the equatorial plane. On the other hand, the types Ib, IIb and IIIb are characteristic for the 'upper' stable halo orbits which occur near the rotation axis. Behaviour of the type IIIc represents an intermediate situation.

Kerr black hole in uniform magnetic field

In order to incorporate the large scale magnetic field in our considerations, we employ Wald's test-field solution [14]. The electromagnetic field is given in terms of the vector potential

A t = 1 2 B 0 g tφ + 2a g tt -1 2 Q g tt -1 2 Q, (38)
A φ = 1 2 B 0 g φφ + 2a g tφ -1 2 Q g tφ , (39)
in the background of Kerr metric

ds 2 = - ∆ Σ dt -a sin θ dφ 2 + sin 2 θ Σ r 2 + a 2 dφ -a dt 2 + Σ ∆ dr 2 + Σ dθ 2 , (40)
where ∆ = r 2 -2r + a 2 , Σ = r 2 + a 2 sin 2 θ (a is the rotational parameter). Q stands for the test charge of the black hole; the terms containing Q may be identified with the components of vector potential of the Kerr-Newman solution [START_REF] Misner | Gravitation[END_REF]. Further, the asymptotic behaviour of components ( 38)-( 39) justifies the identification of the parameter B 0 with the strength of the uniform magnetic field into which the Kerr black hole has been immersed. Wald [14] has shown that in the case of parallel orientation of the spin and the magnetic field B 0 , the black hole selectively accretes positive charges (negative for the antiparallel orientation) until it is charged to the equilibrium value

Q W = 2B 0 a. ( 41
)
We will adopt Q W as a preferred value of the charge in forthcoming examples.

= 1 + (1 -a 2 ) 1/2 . Further classification is given in table 2.
Table 2. Signs of quantities characterizing the motion along stable halo orbits. Values of r and a are taken from the specific regions corresponding to figure 4.

Reg. v h,I Q h,I L h,I A - + + A + - -

Existence and orientation of halo orbits

Our analysis was described in general in section 2. Its application to the case of Wald's test-field is shown in figure 4. Again, we take the advantage of the fact that q appears only in the product with B 0 , which leads us to introduce the effective specific charge Q = qB 0 . Then, the only parameter p ≡ a remains in the relations for v h,i (26), L h,i (28) and Q h,i being expressed by the same formula as q h,i (27). Note that due to the use of Wald's charge, the coefficient A in equation ( 21) vanishes. Moreover, having now the quadratic equations, one of the two possible roots does not take the real values from the interval (-1, 1) in the considered black hole region above the outer event horizon. Thus, there is only one possible root v h,I (r, θ; a) and the related specific angular momentum L h,I (r, θ; a) and effective specific charge Q h,I (r, θ; a) involved in the classification routine. In the region above the outer event horizon, the stable halo orbits are possible only for particles which counter-rotate with respect to LNRF (v h a < 0) and charge of which satisfies the condition Q h a > 0 (regions A, A ). The parallel orientation of the spin a and the magnetic field B 0 (aB 0 > 0) allows counter-rotation along stable halo orbits only for positively charged particles. Moreover, in the potential minimum, there is always v h L h < 0.

Effective potential and classification of halo motion

We can conclude from formulae (40), ( 5), ( 38) and ( 39) that the potential V eff allows several combinations of simultaneous sign reversals which follow from its symmetries. These have to satisfy the conditions sgn(Qa) = const, sgn(aL) = const.

(

) 42
The effective potential (5) does not reveal as many possibilities of stable halo motion as we have seen in the previous case of a magnetic star in section 3.2. Moreover, here, only the topology of the potential is crucial for the classification of the halo motion.

Starting at energy of the potential local minima and rising its level (keeping other parameters constant), we still observe the growing halo lobes in which particles are trapped. The system can be classified in the following way (see figure 5). In the case of the parallel orientation of the spin and the magnetic field, the behaviour of all the presented types determines the motion of positively charged particles counter-rotating along the stable halo orbit. The behaviour of type Ia is characteristic for potentials determining the 'lower' stable halo orbits. Type Ib is characteristic for the potential determining the 'upper' stable halo orbits. Type Ic exhibits an intermediate type.

In all the presented results here, we assume Q = 2B 0 a. But we checked that the stable halo orbits are also possible in the uncharged case, Q = 0, finding no substantial change in the topology of effective potential compared to the charged case.

Discussion: The astrophysical relevance of halo orbits

In this paper, we have concentrated ourselves on aspects of the off-equatorial motion near magnetic compact stars and magnetized black holes. Namely, we discussed the existence of halo orbits, classified them into several categories and investigated their dependence on the parameters. Despite the obvious limitations of our approach, one can recognize several astrophysically relevant features of the halo orbits.

We can list the main limitations of the present approach: (i) Test-particle approximation was employed. This assumes that the medium is highly diluted and the mean free path is comparable with the characteristic size of the system (gravitational radius of the black hole). (ii) Schwarzschild metric was adopted for the external gravitational field of the magnetic compact star and Kerr metric for the rotating black hole, ignoring all other terms that could contribute to the gravitational field. Such contributions could in principle originate from the internal structure of the star and its rotation, or they could be caused by accreted matter outside the central body. (iii) The rotating dipole or the asymptotically uniform structure has been imposed as two examples of the magnetic field and the associated electric field. It was also assumed that the electromagnetic field does not affect the spacetime metric, so it could be treated in the test-field approximation.

The assumptions (ii) and (iii) are less critical because they are well fulfilled under the typical astrophysical conditions. We considered two different examples of the magnetic and gravitational fields, we found that the halo orbits are possible in both. The halo orbits are determined by their location in the poloidal (r,θ) plane, and by values of other parameters, such as the magnetic field strength and the specific charge of the particles. The assumption (i) is the crucial one and it imposes the main restriction on the choice of the system which can be described.

As for the numerical values of the mass-scaled (dimensionless) parameters Ω, a and Q, we can check the adequate choice by looking at the domain of the non-scaled quantities Ω * = Ω/M * , a * = aM * , B * 0 = B 0 /M * , q * = q, R * = RM * and M * , characterizing real objects.

Magnetic star

In the following analysis, we assume masses and radii of magnetic compact stars in the ranges 1 M ≤M * ≤ 3 M and 3M * ≤ R * ≤ 10M * [34], and magnetic fields and rotation frequencies reaching B * 0 = 10 10 T and f * = 10 3 Hz (Ω * . = 6283 rad/s) [35][36][37]. + The limits of the used scaled parameters can be determined from the figure 6, which can be read in the following way. In the upper left plot, starting from the assumed astrophysically relevant values of angular velocity Ω * ≤ 10 3.8 , we can determine values of Ω for a fixed mass M * . For example, light magnetic compact stars with M * = 1 M can rotate with the scaled angular velocity up to Ω . = 0.031, while the heavy ones with M * = 3 M up to Ω . = 0.093. In the upper right plot, within the astrophysically relevant values of magnetic field strength B 0 ≤ 10 10 , we can determine values of B 0 for a fixed mass M * . * We notice that light magnetic compact stars with M * = 1 M can be endowed with magnetic fields of the scaled strength up to B 0 . = 10 -5.4 , while the heavy ones with M * = 3 M up to B 0 . = 10 -4.9 . In the lower left plot, for the scaled magnetic field strength B 0 ≤ 10 -4.9 , we can determine values of M for a fixed scaled radius R. We can also see that small magnetic compact stars with R = 3 allow scaled magnetic dipole momenta up to M . = 10 -3.9 , while the large ones with R = 10 up to M . = 10 -2.0 . In the lower right plot, for the specific charges q ≤ 10 21 , we can determine values of Q for a fixed magnetic dipole moment M. Note that small scaled magnetic dipole momenta of M = 10 -15 allow effective specific charges up to Q = 10 6 , while the big ones of M = 10 -2 can support values up to Q = 10 19 .

+ The precise determination of the mentioned limits is still of high interest, being very complex, dependent on the considered equation of state, etc. Moreover, not all the theoretical calculated models have their observable counterparts. For instance, as for the masses of neutron stars, the narrower range 1.4 M M * 2.5 M is widely stated for the astrophysically relevant neutron stars. As for the radii, most of the realistic equations of state imply the lower limit only R * ≈ 3.5M * [START_REF] Glendenning | Compact stars[END_REF]. On the other hand, the existence of extremely compact stars with R * < 3M * is also discussed. Some models for the so-called Q-stars allow the lower limit even R * ≈ 2.8M * [START_REF] Bahcall | [END_REF][40][41]. In this paper, the limits are considered only for order estimates of the used dimensionless parameters, thus their exact values are of less importance here.

* The value of Ω * is in physical (SI) units, i.e. rad/s. Similarly, B * 0 is expressed in units of Tesla. ). For particles with the specific charge q = 10 12 (charged dust grains), we have the effective specific charge Q . = 21 (see figure 6).

Kerr black hole in uniform magnetic field

While masses of astrophysical black holes are known with relatively good precision and customarily categorized in three groups of stellar-mass black holes (M * 30 M ), intermediate mass black holes (10 2 M M * 10 5 M ), and supermassive black holes (M * 10 6 M ), much less is known about the intensity of cosmic magnetic fields surrounding the black holes. We assume masses of the black holes in the interval 3 M ≤M * ≤ 10 10 M [42,[START_REF] Czerny | Multifrequency Behaviour of High Energy Cosmic Sources[END_REF] and the galactic magnetic field strength up to B * 0 = 10 -6 T. Near supermassive black holes with jets the magnetic field is thought to be still significantly stronger [START_REF] Reynolds | [END_REF].

The limits of the used scaled parameters can be determined from the figure 7, which can be read in the following way. In the left plot, within the assumed astrophysically relevant values of the magnetic field strength B * 0 ≤ 10 -6 , we can determine values of B 0 for a fixed mass M * . We can see that light black holes with M * = 3 M allow values of the scaled magnetic field strength up to B 0 . = 10 -20.9 , while the heavy ones with M * = 10 10 M up to B 0 . = 10 -11.4 . In the lower right plot, within the astrophysically relevant values of the specific charge q ≤ 10 21 , we can determine values of Q for a fixed scaled magnetic dipole strength B 0 . We can see that weak magnetic field strengths of B 0 = 10 -25 allow values of the effective charge up to Q = 10 -4 , while the strong ones of B 0 = 10 -11.4 up to Q = 10 9.6 .

For instance, a rotating black hole with the mass M * = 10 M . = 14750 m in the asymptotically uniform magnetic field of the strength B * 0 = 10 -8 T . = 2.87 × 10 -27 m -1 is described by the parameter B 0 . = 4.24 × 10 -23 . For particles with the specific charge of magnitude q = 10 21 (electrons), we have the scaled effective charge Q . = 0.042 (see figure 7). Realistic values for the scaled spin a are limited by the value a = 1. For the greater values of the spin, the black-hole horizons disappear and the Kerr metric describes the naked-singularity spacetime.

Possible applications

Two kinds of charged particles seem to be relevant with respect to the halo orbits. Firstly, in the case of electrons and protons, the magnitude of specific charges has a large value, and so the halo orbits can exist only very close to the black hole horizon or near an ultra-compact magnetic star. This kind of motion has been recently explored in the context of magnetospherical motion in stellar-mass black hole binaries [45]. Magnetically driven oscillations of matter near black holes have recently been explored also in the context of equatorial motion [33].

Another promising application of the halo motion concerns charged dust grains which typically acquire the magnitude of specific charge much smaller than that of elementary particles. Therefore, the gravitational effects influence the motion of dust grains more efficiently, and so they can form the halo orbits at larger distances from the star surface. Infrared observations are suitable to test the possible existence of halo orbits of dust grains around magnetic compact stars [25]. The debris discs of such particles could perturb the electromagnetic signal of the compact star and give rise to quasi-periodic modulation.

It is also well known that the frequencies of epicyclic geodesic motion [46,47] have a substantial role in the resonant (or other) models of high-frequency quasiperiodic oscillation observed in the black hole [48][49][50] and neutron star [51] low mass X-ray binaries, and their magnetic-field induced modification could be of high importance, as shown in the case of equatorial motion [33].

Conclusions

Our calculations reveal the properties of halo motion of charged particles near magnetic compact stars and black holes. These can be relevant e.g. for charged dust grains when they acquire a small magnitude of electric charge and occur in the strong gravitational field. Such particles can originate from supernova material which falls back and migrates into the pulsar vicinity before evaporating.

The halo orbits can take place in axially symmetric systems: near a magnetic compact star, which we modelled by the Schwarzschild geometry and a test rotating dipole magnetic field, as well as near a Kerr black hole immersed in an asymptotically uniform magnetic field.

In this paper, we concentrated ourselves on the methodological aspects of the halo orbits and their systematic classification. Our main interest concerned the properties of the halo orbits in situations when the motion takes place in a strong gravitational field. We only briefly touched the astrophysical applications and will investigate them further in the following work. Table 3 summarizes the basic characteristics of different examples of strongly gravitating systems which exhibit the halo orbits. Table 3. Orientation of motion of electrically charged particles along stable halo orbits and the corresponding sign of the charge. The results correspond to the configuration where the magnetic field (characterized by the magnetic dipole moment M or magnetic field strength in the equatorial plane B 0 ) is oriented in the same direction as the rotation vector, characterized by the signs of Ω or a.

Compact object

Charge of particles Orientation of motion As for the magnetic compact stars, results of our discussion of the stable halo orbits orientation are qualitatively consistent with those obtained previously in Newtonian and the pseudo-Newtonian investigations [START_REF] Kovář | Classical and Quantum[END_REF]18]. We notice that the strong gravitational field and the presence of event horizon affect the form of the effective potential and introduce significant changes with respect to the Newtonian case. However, we also notice that the pseudo-Newtonian approach can give surprisingly good information on the loci of halo orbits in spherically symmetric gravitational fields of magnetic stars, in which case the results qualitatively agree with the exact general relativistic approach presented here. On the other hand, the pseudo-Newtonian modelling of Kerr spacetimes [52,53] is much more complicated, and seems to be less fruitful for modelling the phenomena in the Kerr geometry, which itself is relatively simple.

The effective potential of the halo motion exhibits a number of qualitatively different types. Two basic categories can be distinguished. The first one allows particles, bounded within the halo potential lobes, to commute across the equatorial plane. On the other hand, the second category does not allow such an interconnection, as the two lobes are entirely disjoint. The off-equatorial halo motion is an interesting phenomenon not only because it has not yet been explored in full detail in the literature, but also because it shows a variety of orbits, depending on orbital parameters: off-equatorial stable circular motion which does not cross the equatorial plane; small oscillations in the radial and vertical directions, with characteristic frequencies around the circular orbits; and the motion that traverses between the lobes when the energy is increased to sufficiently high levels. Recently, the off-equatorial motion of charged particles near Kerr black hole in asymptoticaly uniform magnetic field has been discussed from a different point of view in the ref. [54] as well.

Figure 1 .

 1 Figure1. The range and orientation of particle motion along stable halo orbits near a magnetic star. The star rotates with angular velocity Ω. This figure has been constructed for the latitude θ = π/3 and it represents a typical result for the effective potential V eff (r, θ; L, Q, Ω), based on the behaviour of three orbital velocities v h,I (r, θ; Ω), v h,II (r, θ; Ω) and v h,III (r, θ; Ω) (roots of equation (21)). These correspond to the specific angular momenta L h,I (r, θ; Ω), L h,II (r, θ; Ω), L h,III (r, θ; Ω) and the effective specific charges Q h,I (r, θ; Ω), Q h,II (r, θ; Ω), Q h,III (r, θ; Ω). The gray areas correspond to the parameters for which the effective potentials V eff (r, θ; L h,I , Q h,I , Ω) (left panel), V eff (r, θ; L h,II , Q h,II , Ω) (middle panel) and V eff (r, θ; L h,III , Q h,III , Ω) (right panel) develop local minima. Only the regions under the light cylinder are considered. The thick vertical lines denote position of the event horizon at r = 2. Description of the regions is summarized in table 1.

Figure 2 .Figure 2 .

 22 Figure 2. Behaviour of effective potential V eff for charged particles moving near magnetic star. Qualitatively different types of behaviour of the effective potential occur, with two local minima in the poloidal plane. The lobes are limited by the radius of the star and the light cylinder (thick curves limiting the gray areas). The equipotential surfaces are shown by solid curves, whereas the dot-and-dashed and dashed lines correspond to zeroes of ∂ r V eff and ∂ θ V eff , respectively. The presented types correspond to the following combinations of parameters: Ia: Q . = -5.72, L . = 0.876, Ω . = 0.0115, R = 4 and θ h = π/3, r h = 5 (parameters taken from region A in figure 1); Ib: Q . = -14.6, L . = 0.343, Ω = 0.03, R = 4 and θ h = π/9, r h = 5; IIa: Q . = 5.77, L . = -2.30, Ω = 0.0505, R = 4 and θ h = π/3, r h = 12 (parameters taken from region B in figure 1); IIb: Q . = 25.2, L . = -0.788, Ω = 0.0505, R = 5 and θ h = π/9, r h = 12.

Figure 2 .

 2 Figure 2. (Continued.) The presented types correspond to the following parameters: Ia : Q . = -5.72, L . = 0.876, Ω . = 0.0115, R = 4.85 and θ h = π/3, r h = 5 (parameters taken from region A in figure 1); IIb : Q . = 17.1, L . = -0.617, Ω = 0.0505, R = 5.1 and θ h = π/9, r h = 8.

Figure 3 .

 3 Figure 3. Three-dimensional trajectories of charged particles (upper panels) with Q .= -5.72 and the corresponding poloidal projection of the motion (lower panels), moving in the separated and merged halo lobes in the vicinity of a magnetic star with R = 4 and Ω . = 0.0115. The trajectories were calculated numerically by using the Lorentz equation of motion(8) with the initial conditions L . = 0.876, φ(0) = 0, r(0) = 5, θ(0) = π/3, and E = 0.8481, u θ (0) = 7 × 10 -4 , u r (0) . = 0.021 (motion in separated halo lobe), and E = 0.8485, u θ (0) = 7 × 10 -3 , u r (0) . = 0.020 (motion in merged lobes).

Figure 4 .

 4 Figure 4. The range and orientation of particle motion along stable halo orbits near a Kerr black hole (with the spin a) embedded in a uniform magnetic field. This figure has been constructed for the latitude θ = π/3 and it represents a typical result for the effective potential V eff (r, θ; L, Q, a) based on the investigation of the orbital velocity v h,I (r, θ; a) (root of equation 21 with A = 0). The corresponding specific angular momenta and charges are L h,I (r, θ; a) and Q h,I (r, θ; a). The gray regions correspond to the the local minima of the effective potential V eff (r, θ; L h,I , Q h,I , a). The thick solid curve shows the position of the event horizon at r = 1 + (1 -a 2 ) 1/2 . Further classification is given in table 2.

Figure 5 .

 5 Figure 5. Behaviour of the effective potential V eff for charged particles moving near Kerr black hole in the uniform magnetic field. These are typical examples of behaviour of the potential with two potential minima, again introduced in terms of the equipotential surfaces. The dashed and dot-and-dashed curves correspond to zeros of ∂ θ V eff and ∂ r V eff , respectively. The gray region corresponds to the interior of black hole, limited by the outer event horizon. The presented types correspond to the following parameters: Ia: Q . = 1.37, L . = 3.76, a = 0.6 and θ h = π/3, r h = 2.6; Ib: Q . = 4.79, L . = 2.76, a = 0.9 and θ h = π/9, r h = 3; Ic: Q . = 1.60, L . = 4.02, a = 0.55 and θ h = π/3, r h = 2.5. The parameters of all the presented types are taken from the region A in figure 4.

Figure 6 .

 6 Figure 6. Relations between the physical parameters (magnetic field strength B * 0and angular velocity of rotation Ω * ) and the corresponding scaled (dimensionless) parameters, characterizing a magnetic star and the charged particles moving around it. The dashed lines denote the assumed limits and some representative values. The gray regions are astrophysically irrelevant, within the assumed limitation.

Figure 7 .

 7 Figure 7. As in the previous figure, but for rotating black holes in the uniform magnetic fields.

¶ Because of the complex spectrum of the different types, we present only the most representative cases in the following list.

Acknowledgments

The authors thank M. Urbanec for useful discussions concerning section 5.1. Institute of Physics and the Astronomical Institute have been operated under the projects MSM 4781305903 and AV 0Z10030501, and further supported by the Centre for Theoretical Astrophysics LC06014 in the Czech Republic. JK, VK and ZS thank the Czech Science Foundation (ref. P209/10/P190, 205/07/0052, 202/09/0772). OK acknowledges the doctoral student program of the Czech Science Foundation (ref. 205/09/H033).


George Drettakis

Eugene Fiume

Alain Fournier

Tightly-Coupled Multiprocessing for a Global Illumination Algorithm

A prevailing trend in computer graphics is the demand for increasingly realistic global illumination models and algorithms. Despite the fact that the computational power of uniprocessors is increasing, it is clear that much greater computational power is required to achieve satisfactory throughput. The obvious next step is to employ parallel processing. The advent of affordable, tightly-coupled multiprocessors makes such an approach widely available for the first time. We propose a tightly-coupled parallel decomposition of FIAT, a global illumination algorithm, based on space subdivision and power balancing, that we have recently developed. This algorithm is somewhat ambitious, and severely strains existing uniprocessor environments. We discuss techniques for reducing memory contention and maximising parallelism. We also present empirical data on the actual performance of our parallel solution. Since the model of parallel computation that we have employed is likely to persist for quite some time, our techniques are applicable to other algorithms based on space subdivision.

Introduction

In order to generate realistic images, global illumination algorithms tend to require vast amounts of main memory and computational resources. The demand for increasing visual realism is unlikely to diminish, and any progress made in improved visual realism simply results in increased expectations. Consequently, the performance of implementations of these algorithms on uniprocessors will continue to be unsatisfactory, despite constant improvements in CPU and memory speeds. Parallel solutions to such problems have often been advocated, but until the advent of affordable multi-processors, this option has not been widely available. These multi-processors are becoming quite popular, and it is evident that a popular mode of large-scale computation will be to have many such machines co-existing on a local-area network. Given this newly emerging computational environment, it is of some use and interest to apply it to intensive applications such as global illumination.

In this paper, we shall concern ourselves with a tightly-coupled parallel implementation of FIAT, a global illumination algorithm that we have recently developed [START_REF] Fournier | FIAT LUX[END_REF]. The time and space requirements of FIAT often exceed the capabilities of large uniprocessors. However, it responds nicely to parallelism. While we have naturally tailored our parallel decomposition to the problem at hand, we feel that a number of more general principles arise from our work: that acceptable speedup requires careful reconsideration of the algorithm. that the architecture of the target environment must be taken into account in the parallel algorithm. that first-order approaches to parallelism such as simple locking and buffering achieve suboptimal and sometimes even subserial speedup. that better serial algorithms can result from parallelism.

We believe that the basic tenets of our approach readily transfer to other intensive graphical applications.

In the next section, we give a brief overview of FIAT. Further details of the algorithm will be presented as needed. We shall discuss a specific parallel implementation of FIAT, and the effect of the implementation on the performance of FIAT. It is unfortunately not common for parallel approaches to progress beyond the paper-design stage. In our case, we are able to report on the performance of a true implementation. Lastly, we draw conclusions and examine further extensions to our work.

An Overview of FIAT

FIAT is a physically-based approach to the simulation of global (or environmental) illumination. Related approaches include radiosity [START_REF] Goral | Modelling The Interaction Of Light Between Diffuse Surfaces[END_REF], ray-tracing [START_REF] Whitted | An Improved Illumination Model for Shaded Display[END_REF], and hybrid two-pass schemes ([Imme86], [START_REF] Shao | A New Radiosity Approach by Procedural Refinements for Realistic Image Synthesis[END_REF], [START_REF] Sillion | A General Two-Pass Method Integrating Specular and Diffuse Reflection[END_REF]). Our approach differs in that it directly handles a wide class of light sources together with surfaces having arbitrary reflectance behaviours. This can be contrasted with traditional radiosity-based approaches which typically are constrained to diffuse reflectors and emitters, or ray-tracing approaches, which focus on specular reflectors and point or parallel light sources. FIAT directly handles the gamut of reflectance functions from the purely diffuse to the purely specular. A detailed presentation of FIAT is given in [START_REF] Fournier | FIAT LUX[END_REF] and also in [START_REF] Ouellette | Modelling Light Propagation and Global Illumination[END_REF]. Our goal in this section is to describe very briefly the overall operation of FIAT, so that its parallel solution can be understood, and so that the application of our techniques to other algorithms is clear. FIAT is still somewhat experimental, but recent results have demonstrated novel illumination phenomena, such as secondary specular reflection and caustics -a phenomenon that is very difficult to achieve using other means (see Figure 1(a)).

FIAT takes as input a description containing geometric and light-source data. An 'initial cube' that bounds the scene is computed, and each face of the cube is tessellated into the same number of side elements or sexels. Each sexel occupies a small square on the side of a cell, and contains light information in the form of parallel beams of square cross-section. Each parallel beam represents a sample of light power emanating from the sexel in a specific direction. Many directions are used to sample a hemisphere of incoming light distribution. The goal of FIAT, coarsely speaking, is to convert global light information into a large set of sexels distributed throughout the environment. The sexels in turn contain many samples of light flux in various directions (beams). The initial cell, the sexels and the beams are all shown in Figure 2.

The reduction is achieved by maintaining a power budget for each cell, which takes into account incoming, outgoing, emitted, and absorbed light power. A cell is balanced if the incoming and emitted light power totals the outgoing and absorbed power. There is initially just one cell, the bounding volume, and the process begins by computing the power entering the volume from external light sources. This results in a net power imbalance, which FIAT will try to rectify. If a cell is empty or contains objects that can be dealt with directly, then a new power balance is determined. If the cell cannot be dealt with directly, then it is subdivided in octree fashion, and relevant sexel information is transferred to the subcells.

We consider in slightly greater detail the case in which a cell can be dealt with. If a cell is empty, we simply propagate every incoming light beam to relevant outgoing sexel. If the cell is not empty, an intersection computation between a specific light beam and the scene within the cell is performed. Light beams that do not intersect with any objects within the cell are propagated as in the empty-cell case. For those beams that do intersect with an object, a new power distribution about the point of intersection is computed by invoking a local illumination model, with the beam acting as a light source. The outgoing distribution may be quite complex, since it is not assumed that the reflectance function of the intersected object is perfectly diffuse or specular. This distribution is sampled in many directions, and the relevant outgoing sexels are updated. If the intersection point corresponds to a point that is visible from the eye position, then the shade of this point in the direction of the viewpoint is computed and added to a relevant pixel entry in a light update buffer or LUB. Thus an image is progressively produced as a function of cell

One Subdivision

Detail of Sexels on Cell Faces Sexels on a Face Initial Cell Beams Over One Sexel

The Goals for a Tightly-Coupled Parallel Solution

The task of FIAT inside a cell is complex and can involve a large amount of computation and data. To take maximum advantage of the capabilities of parallel hardware we set the following goals: (a) To perform as much of the computation as possible in parallel.

(b) To distribute the computation evenly among all processors.

(c) To find a solution that minimises the amount of locking necessary, and allows the processors to compute independently as much as possible.

(d) To introduce a solution that can be scaled to a large number of processors. The range of processors in which we are interested is that of existing common-bus shared-memory architectures, incorporating from 2 to 16 processors and not more than about 30. For larger-scale parallelism, we prefer to take a hierarchical approach in which several tightly-coupled processors are employed. This approach is an orthogonal dimension to the techniques described in this paper [START_REF] Drettakis | Hierarchical Parallelism for a Global Illumination Algorithm[END_REF].

To achieve these goals we shall first analyse the serial algorithm to determine the computationally expensive sections. We investigate the structure of the algorithm and the data dependencies involved, to decide how much computation in these sections can be parallelised, thus achieving goal (a). Goal (b) is achieved by investigating the ways in which the computation can be segmented into independent components so that the work is evenly distributed. To avoid locking it may be necessary to use temporary memory private to each processor, and examine whether such a solution is actually preferable. For both (b) and (c) in combination with private buffering schemes, changing the order of computation proves to be necessary. Finally the selection of the scheme will be such that goal (d) is satisfied by design.

In the remainder of the paper we present the serial algorithm in more detail, outline a set of test cases which we use to analyse the behaviour of the algorithm, and investigate the choices of parallel schemes.

We present one solution that satisfies the goals to the largest extent, and discuss some of the problems that were dealt with in the design and implementation of this solution. Timing results on three test scenes are presented showing the effect of varying several parameters affecting performance. Some possible alternative schemes for parallelism are discussed, and we close with a summary and some possible enhancements.

Analysis of the Serial Algorithm

The computation in a cell can be split into two components. The first, called balancing, involves distributing light entering a cell as a result of interaction with objects in the cell; also involved is the propagation of light beams within the sexels on one face of a cell to the sexels on some other face if the beam does not intersect with anything. The second component, shading, involves updating the section of the LUB associated with the cell. This is done by determining the visible surface for each pixel and shading based on the incoming light that impinges on this visible point. These two processes are independent, but use some common subroutines.

The serial algorithm.

For the purposes of this discussion we present a high-level pseudo-code version of the balancing (Figure 3

Analysing the behaviour of FIAT in a cell.

From the pseudo-code presented in Figure 3 it is evident that a large proportion of the computation will be spent in the functions that are executed in the inner loops. For balancing, these are the beam intersections and, if distribution is needed, the power calculations and the code that finds the destination sexel. For shading, the corresponding expensive operations are finding the sexels that shine on the visible point, and performing power calculations to determine the shade. To confirm these intuitive estimations profiling was performed for test cases. These test cases will also be used to evaluate the performance of the parallel solution proposed.

Test Cases

FIAT is extremely flexible in the range of light sources and surface properties that can be simulated. Test cases must therefore be chosen with care to exercise as many modalities as possible. We have employed the following three test scenes.

Scene 1: A shiny ashtray (60% specular, 40% diffuse) lit by a directional light source that is at a 45 degree angle. The model consists of 600 triangles, and as a result many cells contain more than one object. We employ directional light sources to illuminate this environment, which results in the creation of caustics from concentration of light due to reflection. Figure 1(a) shows this scene computed at higher subdivision to demonstrate the effect more clearly.

Scene 2: A room that has 3 walls, a red left wall, a blue right wall, a grey back wall, and a grey ceiling and floor. There are 12 triangles in this scene. All surfaces in this scenes are mostly diffuse, and absorb some of light falling on them. The scene is lit by a extended area light source that is where a front wall would be if the room was closed all around. This scene shows colour bleeding, similar to that shown by radiosity on the ceiling from the two side walls. Figure 1(b) shows this scene computed using the reconstruction techniques that have been developed for FIAT, eliminating some of the artifacts.

Scene 3: The same as Scene 2, with the difference that the floor is now shiny (60% specular 40% diffuse) and acts as a mirror. The image is similar to that presented in Figure 1(b), with a different floor.

For scenes 2 and 3 only a small number of sexels and cell subdivisions are employed. Because memory requirements for these scenes can exceed real memory on some machines, it is occasionally necessary for the purposes of analysis to use scene subsets. For scenes 2 and 3 the subset scenes have the floor and the back wall removed, and for scene 3 the ceiling is shiny instead of the floor that is no longer present. We refer to these scenes as scenes 2a and 3a respectively.

Profiling results.

By using standard profiling tools we have determined which functions for each corresponding part of the balancing and shading computations are the most expensive, and are therefore prime candidates for parallelism. The profiling results confirm what was expected from examining the structure of the algorithm. These results show that overall more time is spent balancing than shading. In all cases 50-60% of the execution time is spent in BalanceCell. Of this time, for Scene 1, about 40% is spent in applying the local illumination model at the intersection point to determine the power outgoing in each direction ( PowerCalc and power assignments). The additional 10% is spent determining which sexel is affected. For Scenes 2 and 3 the time is split almost equally between power calculations and determining the destination sexel. This discrepancy is due to the nature of the ashtray surface, that is highly specular, thus reflecting light in a limited number of directions. The beam intersection operation is also somewhat expensive (2-3% of the computation time).

For shading the time is almost equally split between finding the sexel that affects the visible point being considered ( FindHitSexel) and the time spent applying the local illumination model to shade the point ( ShineBeam). This is true for all scenes.

Linearity of Light Interaction

In order to allow progressive updates to the light update buffer (LUB) that do not depend on the order of application, we have assumed that light-object interaction is a linear operation [START_REF] Fournier | FIAT LUX[END_REF]. For most materials and light sources, this is a perfectly valid assumption. A consequence of this assumption is that light beams can be propagated or distributed in any order, yielding the same cumulative result. This is very important since it facilitates the reorganisation of a large amount of the computation in ways that are more convenient for parallelism.

A Parallel Decomposition Strategy

To devise a good parallel decomposition strategy, one must identify (potentially) non-interfering portions of the computation. It is often the case that although independent tasks exist, parallel execution can result in corruption of some data structures that tasks access at the same time. The obvious first-order approach is to require that a task acquire a lock before entering the critical region of code that modifies common data. The gain in parallelism in this case is often suboptimal. Another solution involving more intensive analysis is to reorganise the computation such that tasks no longer write the same memory locations, either by using a private buffer in each processor or by eliminating the need for such a write. Finally a scheme must be developed that actually assigns these parallel tasks to processors in a manner that keeps each processor as busy as possible at all times. In discussing possible parallel strategies we will address balancing and shading since they are independent computations.

Parallel Balancing.

For each incoming beam in a sexel we have to perform an intersection with the objects in the cell. The intersection operation of an incoming beam is independent of any other intersection of the same form, since no data is written apart from a temporary structure that contains the hit point information if such a point exists. Thus, these intersections can be easily performed in parallel.

Most of the operations involved in the distribution of a light beam about an intersection point are independent. These operations include invoking the local illumination model, finding the sexels to update, and power computations. The complicating issue is that more than one intersection point may contribute to the same outgoing direction of the same sexel (see Figure 4), since there is only one representative of each direction for each sexel.

A more subtle problem appears if we choose to parallelise the actual distribution of light around an intersection point. In DistribLight (Figure 2(a)) we see that at every intersection point a power calculation occurs for each discrete direction in the sphere of directions placed on that point. An obvious way to parallelise is to perform the power calculations for each direction in parallel. At a first glance this seems a reasonable strategy. However FIAT requires that the local illumination model be normalised such that the total power leaving an intersection point is no larger than the power entering it. In other words if the power impinging at point P from a given direction was P imp , and the model used assigns P out i power to be outgoing in direction i (Figure 5), we wish the total outgoing power P out (if d is the number of directions in the sphere around P ):

P out = i =0 Σ d P out i ≤ P imp .
In any physically-based local illumination model, power must be conserved. Unfortunately, few models in computer graphics have this property, and so a pre-normalisation step is required. In the initial serial implementation of FIAT, each P out i is multiplied by a scaling factor S defined as follows. (1)

Consequently,

P out = i =0 Σ d S P out i = S i =0 Σ d P out i = i =0 Σ d P out i P imp i =0 Σ d P out i = P imp .
Normalisation guarantees conservation of power even for unrealistic illumination models. The final value assigned to the destination sexel is S P out i .

In the serial case, the value S can be computed on-the-fly, since one processor handles all the P out i . This is clearly not the case in a parallel solution, so one must arrange to precompute a representative S . Otherwise some processors may have to wait for others to complete before being able to apply S to their outgoing power values. Later, we shall see how a good value for S can be precomputed, so that this potential barrier is removed entirely. We stress that this problem only arises when local illumination models that do not conserve power are employed.

First-Order Parallel Balancing

Without changing the structure of the initial algorithm, we can statically divide the processing of faces, sexels or beams among processors. In such a scheme the main loop of BalanceCell() is executed on every processor but only certain of the faces, sexels or beams are dealt with by each. Parallelism at the face or sexel level is unsuitable for two main reasons. First, such an approach does not scale well (there are 6 faces, and in the worst case there are 6 sexels in a cell). Second, it does not provide very good load balancing as there may be one face or only a small number of sexels that contain all the incoming light. At the beam (direction) level, things are much better as there are usually many beams (a typical instance of FIAT has 192 such beams per sexel), and they are usually evenly distributed. However such a scheme may suffer from the contention problem described in the previous section and illustrated in Figure 4. We briefly describe three ways of addressing this problem without significantly altering the overall structure of the initial serial algorithm.

Locks. Locking critical sections is straight-forward, but can cause excessive waiting. Moreover, the use of many locks (as would be the case in FIAT if each sexel were to have a lock) would pose high additional memory requirements.

Private buffers: Perform all sexel updates in a private buffer, acquire a lock when the buffer is full, and write the updates to the real sexels. This also suffers from the problem of potentially long delay on locks as the lock must be acquired for the full duration of the buffer write-back. A severe practical drawback is that in a memory-intensive application like FIAT, real sexel information migrates to distant locations in memory while buffers are filling. In systems that use look-aside buffers for memory accesses, processor faults are caused at an excessive rate every time the buffer has to be written back (typically thousands of times a second). When this problem appears, the actual performance of this approach can even be subserial.

Pooled Buffers. A third solution is to employ a two-pass approach. Each sexel is preassigned to a processor in some static manner. If there are n processors, each processor maintains n buffers. In the first pass each processor propagates and distributes light, and when a sexel must be updated, it inserts the update information in the buffer corresponding to the processor assigned to this sexel. In the second pass each processor runs through the buffers containing updates for the sexel assigned to it on all the processors, performing the real updates. This solution avoids the problem of frequently accessing distant memory as the set of sexels that each processor accesses is smaller; during the first pass only local accesses are required. This solution yields relatively good results but is limited by the buffer-management and synchronisation overhead.

From the above discussion we see that simple solutions that do not alter the original structure of the serial algorithm or have marginal changes are limited due to synchronisation, memory access or locking constraints. Shortly, we shall describe a solution that allows these problems to be circumvented so that satisfactory speedup is achieved.

Parallel Shading

Shading presents much less difficulty than balancing. The reason for this is that the only data being written in the shading process are the LUB pixels. Therefore if a scheme is found that assigns each pixel to a processor, there is no contention, and each processor can process the pixels assigned to it independently without requiring any synchronisation or locking.

A Fast Low-Contention Parallel Balancing Scheme

To obtain maximal parallelism and minimal memory contention, reorganisation of the computation from the initial serial algorithm is required. Both propagation and distribution of light around an intersection point are based on a global sphere of directions. In the case of propagation, the set of directions within a sexel are defined such that a beam originating from one sexel always updates the beam of a sexel on a different face that has the same direction as the source beam (see Figure 6(a)). In the case of distribution around an intersection point P , the power outgoing in direction i is assigned to a sexel (or sexels) in direction i (Figure 6 Each direction is assigned to a different processor, making this processor responsible for all updates in this direction on any sexel, as a result of either propagation or distribution.

For propagation, this scheme is simple, since the serial counterpart can remain intact with each processor simply working only on those beams that are assigned to it. For distribution of light around intersection points it is necessary for the all intersection points to be available to all processors. In this way each processor can perform the computation for the directions assigned to it for all intersection points.

To allow this, a two-phase buffering scheme is used. In the first phase each processor performs the intersections of the beams it is responsible for, saving the intersections in a private buffer residing in shared memory. When the buffers fill, the processors synchronise and the second phase begins. At this point, each processor runs through the buffers corresponding to all the processors and for every intersection performs the power calculations only for the directions assigned to it and updates these directions on the relevant sexels.

As described earlier, we must ensure that the local illumination model is properly normalised before distribution, for otherwise all processors would have to synchronise at every intersection point to use the total outgoing power for the scaling factor. This would create unnecessary waiting and would severely degrade performance. Fortunately, a good candidate for a normalisation constant is at hand. We shall derive an upper bound for the amount of power that can possibly leave a point P on a given surface if incoming power P imp impinges on P . This upper bound was developed for other purposes in [START_REF] Ouellette | Modelling Light Propagation and Global Illumination[END_REF], but however is exactly what is necessary to overcome this problem. Assuming a Blinn-like reflection model for a surface with constant n > 0, an upper bound P u on outgoing power is:

P u = 0 ∫ 2π 0 ∫ π/2 cos n x dx dy = 2π 0 ∫ π/2 cos n x dx.
We can use P u in place of (2)

one memory servicing all the requests there is an effect of serialisation when memory requirements are acute. This fact is confirmed using the system monitoring facilities that report look-aside buffer faults in the range of 3000-4000 faults/second. This is one of the problems that affects the degradation of efficiency as the number of processors grow.

The following graphs present speedup by the number of processors from 2 to 4. Graph 1 shows the speedup for the scenes 1, 2 and 3 and graph 2 for scenes 2a and 3a. We see from these that the speedup for the small scenes (2a, 3a) is better, and there is little degradation as the number of processors grow. The following two graphs show processor utilisation for the two sets of scenes, again scenes 1, 2 and 3 that have larger memory requirements (Graph 3), and scenes 2a and 3a in Graph 4. From these graphs, which use the same scale for purpose of comparison, observe that the penalty for each extra processor for the large memory cases is 6-12% while for the scenes with low memory requirements it is much less varying from 2-5%. This means that if the memory requirements are smaller the solution can scale to a larger number of processors.

Conclusions and Future Directions

We have presented an approach that allows the use of parallelism for FIAT, a complex global illumination space-subdivision algorithm, and have successfully demonstrated significant speed-up with the addition of processors. The techniques readily migrate to other parallel environments: we have recently migrated the parallel implementation of FIAT to a 7-processor DEC Firefly. We also believe that our techniques to reduce memory contention and locking can be employed in other applications.

Proceedings of Eurographics 90

Montreux, Switzerland, Sept. 1990

FIAT has evolved since the approach described in this paper was originally developed. Reconstruction techniques are used in both the propagation and distribution phase and in the shading, to allow a better representation of light flux using the discrete samples. The methods described in this paper are directly applicable to the new technique, with one exception: there is no longer a one-to-one correspondence between the originating sample (direction) and the destination sample, since more than one sexel and more than one direction are affected. This problem can be overcome with some additional effort.

In [START_REF] Drettakis | Hierarchical Parallelism for a Global Illumination Algorithm[END_REF] a full hierarchical solution is designed and implemented. The octree structure used to represent the cell subdivision is distributed across a network of parallel processor stations. The tightlycoupled solution described here is run on each such station, and special algorithms and techniques are developed to allow maximal benefit by distributing both memory and computation.

Figure 2 :

 2 Figure 2: Initial Cell, Subdivision, Sexels, Beams balancing. The algorithm continues until the overall power imbalance falls below a threshold value.

  (a)) and shading algorithms (Figure 3(b)) for the serial version of FIAT, shown on the following page.

  Figure 3 (a) The serial algorithm for balancing.

Figure 3

 3 Figure 3(b) The serial algorithm for shading.

Figure 5 :

 5 Figure 5: Outgoing power in direction i

Figure 6 :

 6 Figure 6: Assigning directions.

  in equation (1), giving us a revised value of S :S = P u P imp .

  . . . . . . . . . . . . . . . . . . . . . . . . . . Scene 3a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Graphs 1 and 2: Speedup for Scenes 1, 2 and 3 and for Scenes 2a and 3a.

  . . . . . . . . . . . . . . . .. .. .. .. . .. .. .. .. ... Scene 3a . . .. . . . .. . . . .. . . . . .. . . . .. . . . .. . . Graphs 3 and 4: Processor utilisation for Scenes 1, 2 and 3 and for Scenes 2a, 3a.

Figures 1

 1 Figures 1(a)-(b). Figure 1(a) depicts Scene 1 calculated at high subdivision, showing a caustic in the centre of the ashtray. Figure 1(b) depicts Scene 2, calculated using the new reconstruction techniques for both distribution and shading. Colour bleeding is clearly shown.

Different Intersection Points Update Same Beam.

	P 1	. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . P 2	sexel boundaries
	Figure 4: S =	i =0 Σ P out i d P imp .

† The authors gratefully acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada, the Information Technology Research Centre of Ontario, Apple Computer Inc., and the University of Toronto.

P u can be precomputed once for each distinct n using a numerical integration package, or an analytical derivation if n is integral, and stored with the surface description. To perform normalisation each process needs only to apply equation (2), and synchronisation is avoided.

We summarise our approach by the following pseudo-code. On each processor the code shown in Figure 7

Goals Achieved by the Parallel Approach

The solution just described achieves the goals set out at the beginning of this paper:

All expensive portions of the computation are performed in parallel (i.e. beam intersections, power calculations, destination sexel determination for distributing).

The scheme has a satisfactory 'natural load balancing' property, because load distribution is tied to distribution of directions.

Memory contention is minimised and locks are almost entirely avoided. This solution also has the benefit that there are very few barriers at which the processors have to synchronise. The only required synchronisation point occurs when the intersection buffer fills.

The solution can in principle be scaled to the number of processors equal to the number of directions. However the benefit of this order of parallelism would quickly be overcome by the overhead of maintaining and synchronising a large number of processors. Although we report results on machines consisting of 4 processors or fewer, it seems that any number of processors that fit on a single-bus system (less than 32 in most existing systems) would perform well.

From the above we see that from the point of view of design this scheme has the desirable properties that are required from the parallel balancing scheme. As we shall see in the section on timing results the solution performs satisfactorily on a real parallel system.

A Parallel Approach to Shading

From Figure 3(b) we can see that parallelising the shading computation is an easy task. All that is required is to assign each pixel to a different processor that then performs the ray intersection and the power calculations when collecting light from the relevant sexels in parallel. If we have n processors, the parallel shading algorithm executed on processor i in parallel is shown in Figure 8. For every scan line each processor with processor-id i shades pixels i + jn where n is the number of processors and a ≤ j ≤ b such that a = min k : xmin ≤ (i + kn ) and b = max k : (i + kn ) ≤ xmax .

The shading scheme successfully achieves the goals put forward:

Shading expense is determined by light collection, which requires a power calculation for the light shining from the relevant sexels, finding these sexels and in the intersection of the ray from the eye. All are performed in parallel.

The load balancing is dependent on the values of xmin, xmax, ymin, ymax.

This scheme requires no locking or synchronisation for the duration of the shading computation.

The parallel shading approach we have proposed scales theoretically to as many processors as there are pixels in a region of the LUB related to a cell, although diminishing returns would likely set in well before this upper bound is reached.

Effect of the Parallel Solution on the Serial Algorithm

As we tested the parallel implementation, we discovered an interesting result. The parallel algorithm was tested initially for various test cases that fit in the real memory of the systems used. As soon as subdivision levels in FIAT are raised the memory used quickly becomes much larger than real memory, therefore incurring the overhead of paging. Surprisingly, comparisons with the original algorithm appeared to be giving the parallel version speedup 2-2.5 times the number of processors. For example, on a 4 processor system the speedup was 9.8. It quickly became clear that the reorganisation of the computation by buffering intersection points and separating the direction space into sections equal to the number of processors provided better memory-access patterns. By simulating the parallelism in the uniprocessor case the serial version achieved much better performance when paging overwhelms the computation. Thus our parallel algorithm improved serial performance as well.

Timing Results

Tables 1-5 show the timing results for the three selected test scenes. Also shown is speedup and processor utilisation (efficiency), and scaled speedup and scaled processor utilisation. Let t n be the running time of the implementation on n ≥ 1 processors. We define the speedup S n = t n / t 1 , and the processor utilisation PU n = S n / n . The serial time reported in the results is the time in which the parallel algorithm must run in serial on one processor, and includes things such as the time spent reading in the scene, which is obviously large for the 600 triangles in scene 1. The serial time is the average of the reported times spent in serial mode for each run that differed by a maximum of .5 seconds. If s is this serial time we define scaled speedup as S′ n = (t ns ) / (t 1s ) and scaled processor utilisation is PU′ n = S′ n / n . Tables 1-3 were computed on an IRIS 4D/240 with 4 16MHz MIPS R3000 processors and 64 Megabytes of real memory, while the results reported in Tables 4 and5 were computed on an IRIS 4D/120S with 4 10MHz MIPS R2000 processors and 16 Megabytes of main memory. The implementation is in 'C' using the system-provided parallel support library routines for parallel tasks and semaphores. The speedup and efficiency achieved are generally satisfactory. There are two factors affecting speedup that can be noted from the results.

Light source and surface types. Our solution depends on the type of light sources and the type of surfaces existing in the scene. The results for scene 1 are worse overall, since the light source is directional and the surfaces are specular. This means that a limited number of directions have more power than others and as a result processor utilisation is suboptimal.

Large memory requirements and hardware. In scenes 1-3 we see worse performance than scenes 2a and 3a. The reason for this is the memory system hardware of the specific systems used. The special fast look-aside buffers are strained with the memory requirements of FIAT in each cell, and since there is only