text:
# Dynamics of Relativistic Flows

## 1 Introduction

A significant feature associated with quasars and active galactic nuclei is the existence of relativistic outflows known as “jets.” The jets frequently come in pairs that originate from a central “engine” that is conjectured to be a massive rotating black hole surrounded by an accretion disk. Rather similar phenomena are also observed in most of the X-ray binary systems in our galaxy; therefore, these are usually called “microquasars.” It is thought that the jets are emitted along the rotation axis of the Kerr black hole. The double-jet structure is consistent with the symmetry of a Kerr source under reflection about its equatorial plane as illustrated in Figure 1. The observational aspects of the motion of astrophysical jets are discussed in .

The purpose of this paper is to discuss a relativistic tidal effect that can be used to predict the speeds of astrophysical jets relative to their ambient media. The theory of the formation and propagation of jets per se is beyond the scope of this paper. We will consider a Kerr black hole with plasma outflow along its rotation axis and study the dynamics of a clump in the flow *after it has emerged from its source environment*. To explore the transient role of the novel relativistic tidal effect near the source, we assume that large-scale electromagnetic fields are dynamically important. We thus ignore the complicated plasma processes that actually determine the electrical drift velocity of the ambient medium. Our approach may be contrasted with the purely hydrodynamic models of relativistic outflows: all force-free as well as MHD wind models assume at the outset that the plasma flows at the drift velocity. It is important to emphasize the idealized nature of the model employed here; our purpose is simply to illustrate the regime in which previously ignored relativistic tidal forces are significant for the dynamics of relativistic outflows.

The mechanism for clump formation and the possible values of the clump’s initial speed as well as other characteristics are beyond the scope of our paper. Once emitted, the motion of the clump relative to the ambient medium is subject to the ubiquitous tidal forces of the central source together with the Lorentz force due to the large-scale electromagnetic field that is assumed to be present in the ambient medium. In addition to various other simplifying assumptions throughout this work, all radiative as well as nonlinear plasma effects  are also neglected here.

With respect to a background Fermi frame that is fixed on the ambient medium with coordinates $`(T,𝐗)`$ such that $`𝐗=(X,Y,Z)`$ and the $`Z`$-axis represents the axis of rotation of the black hole, the equation of motion of a plasma clump along the rotation axis relative to the ambient medium is given by a “Newtonian” equation of motion, involving the superposition of external accelerations, of the form

$$\frac{d^2Z}{dT^2}=A_T+A_L,$$
(1)
where $`A_T`$ is the relativistic tidal acceleration due to the central rotating mass and $`A_L`$ is the Lorentz acceleration due to the electromagnetic field of the ambient medium. The nature of these accelerations are discussed in turn in Sections 2, 3 and 4. The main dynamical features of the resulting equation of motion are described in Section 5. Section 6 contains a brief discussion of our results. We employ units such that $`c=1`$.

## 2 The Generalized Jacobi Equation

Imagine a background global inertial system with coordinates $`(T,𝐗)`$ and two neighboring test particles in free fall in a Newtonian gravitational potential $`\mathrm{\Phi }(𝐗)`$. Let $`𝐗_1(T)`$ and $`𝐗_2(T)`$ denote the paths of the two particles. Neglecting the gravitational attraction between the test particles, their relative position $`𝐗(T)=𝐗_1(T)𝐗_2(T)`$ is determined by the tidal equation

$$\frac{d^2X^i}{dT^2}+K_{ij}(T)X^j=0,$$
(2)
where $`K=(K_{ij})`$ is a symmetric matrix given by

$$K_{ij}=\frac{^2\mathrm{\Phi }}{X^iX^j}$$
(3)
evaluated along the orbit of one of the particles designated as the reference trajectory. Note that $`\mathrm{tr}(K)=^2\mathrm{\Phi }=4\pi G\rho _N`$, where $`\rho _N`$ is the Newtonian mass density of the source. If the motion occurs exterior to the source, then the matrix $`K`$ given by (3) is trace-free and harmonic. The tidal acceleration in (2) has been evaluated to linear order in the separation of the particles.

In the theoretical framework of general relativity, the free neighboring test particles follow timelike geodesics of the spacetime manifold. Their relative motion is described by the geodesic deviation equation. This equation may be expressed in terms of a Fermi normal coordinate system that is constructed along the reference geodesic. A Fermi frame is an almost inertial coordinate system in a finite cylindrical region along the worldline of the reference geodesic . We now let $`(T,𝐗)`$ refer to the Fermi coordinates in such a system, where the spatial origin $`𝐗=0`$ is occupied by the fiducial test particle. Then, the Jacobi equation (2) describes the relative motion of the other test particle in Fermi coordinates and

$$K_{ij}={}_{}{}^{F}R_{0i0j}^{},$$
(4)
provided that the speed of relative motion is negligibly small compared to $`c=1`$. Here

$${}_{}{}^{F}R_{\alpha \beta \gamma \delta }^{}=R_{\mu \nu \rho \sigma }\lambda _{(\alpha )}^\mu \lambda _{(\beta )}^\nu \lambda _{(\gamma )}^\rho \lambda _{(\delta )}^\sigma ,$$
(5)
so that the Riemann tensor in the Fermi frame evaluated along the reference trajectory $`(T,\mathrm{𝟎})`$ is in effect the Riemann tensor projected onto the nonrotating orthonormal tetrad $`\lambda _{(\alpha )}^\mu `$ that is carried along the fiducial geodesic and upon which the Fermi system is constructed. Specifically, the temporal axis of the fiducial observer $`\lambda _{(0)}^\mu =dx^\mu /d\tau `$ is the vector tangent to its worldline, $`\tau `$ is its proper time and $`\lambda _{(i)}^\mu ,i=1,2,3`$, are its spatial axes. One can show explicitly that in the Newtonian approximation of general relativity (4) reduces to (3), so that the formal analogy developed here has a deep physical basis.

The deviation between geodesics is taken into account only to linear order in the Jacobi equation; in fact, higher-order terms can be neglected so long as the deviation is very small compared to the radius of curvature of spacetime $``$. On the other hand, in certain circumstances, the relative speed of neighboring geodesics may not be small compared to unity. The corresponding generalization of the Jacobi equation is given by

$$\begin{array}{c}\frac{d^2X^i}{dT^2}+{}_{}{}^{F}R_{0i0j}^{}X^j+2{}_{}{}^{F}R_{ikj0}^{}V^kX^j\hfill \\ \hfill +(2{}_{}{}^{F}R_{0kj0}^{}V^iV^k+\frac{2}{3}{}_{}{}^{F}R_{ikj\mathrm{}}^{}V^kV^i+\frac{2}{3}{}_{}{}^{F}R_{0kj\mathrm{}}^{}V^iV^kV^{\mathrm{}})X^j=0.\end{array}$$
(6)
This generalized Jacobi equation has been discussed in detail in . If we limit our attention to motion along a fixed direction in space, say along the $`Z`$-direction that characterizes the rotation axis of a Kerr system in the Fermi frame, then (6) reduces to

$$\frac{dV}{dT}+\kappa (12V^2)Z=0,$$
(7)
where $`V=dZ/dT`$ and $`\kappa ={}_{}{}^{F}R_{0303}^{}`$. This equation has an interesting feature that was first pointed out in : Special solutions of (7) exist that correspond to uniform rectilinear motion with limiting velocities of $`V_T=\pm (2)^{\frac{1}{2}}`$. For $`\kappa <0`$, these correspond to attractors if

$$_{T_0}^{\mathrm{}}\kappa (T)(T+C_0)𝑑T=\mathrm{}$$
(8)
for an arbitrary real number $`C_0`$ and an arbitrary initial time $`T_0`$ .

It is useful to consider the application of (7) to motion along the rotation axis of a Kerr source. The main results are expected to hold qualitatively for any central rotating and axisymmetric configuration such as, for instance, the Kerr spacetime endowed with an infinite set of higher moments . The motion of free test particles along the symmetry axis of an exterior Kerr spacetime representing the stationary gravitational field of a source of mass $`M`$ and angular momentum $`Ma`$ is given by

$$\left(\frac{dr}{d\tau }\right)^2=\gamma ^21+\frac{2GMr}{r^2+a^2}$$
(9)
in Boyer-Lindquist coordinates. Here $`\gamma >0`$ is an integration constant; in fact, in the $`\gamma 1`$ case it has the interpretation of the particle’s Lorentz factor at infinity $`(r\mathrm{})`$. We take the reference geodesic to represent the motion of the ambient medium with $`dr/d\tau >0`$ along the $`Z`$-axis in the exterior Kerr spacetime and we construct a Fermi coordinate system $`(T,X,Y,Z)`$ in the neighborhood of this geodesic such that as the reference geodesic is approached, $`(T,𝐗)(\tau ,\mathrm{𝟎})`$. Evaluating the Riemann curvature tensor in this frame, we find that $`\kappa k`$ is given by

$$k=2\frac{GMr(r^23a^2)}{(r^2+a^2)^3}.$$
(10)
That (10) is explicitly independent of $`\gamma `$ is a consequence of the fact that Kerr spacetime is of type $`D`$ in the Petrov classification and hence the axis of rotation provides two special tidal directions for ingoing and outgoing trajectories . Moreover, $`a=0`$ corresponds to radial motion in the exterior Schwarzschild spacetime. For a Kerr black hole $`0<aGM`$, where $`a=GM`$ corresponds to an extreme Kerr black hole. To emphasize that our treatment involves motion along the rotation axis of a Kerr system, it is useful to introduce Cartesian coordinates $`(x,y,z)`$ defined in terms of Boyer-Lindquist coordinates $`(r,\theta ,\phi )`$ as usual: $`x=r\mathrm{sin}\theta \mathrm{cos}\phi `$, $`y=r\mathrm{sin}\theta \mathrm{sin}\phi `$ and $`z=r\mathrm{cos}\theta `$. In Section 5, we will employ this $`z`$-coordinate along the Kerr rotation axis that in the Fermi system becomes the $`Z`$-axis. Moreover, we will assume that the reference particle at $`Z=0`$ is such that $`\gamma =1`$ and $`r>\sqrt{3}a`$ along its geodesic path. The $`\gamma =1`$ condition implies that the reference particle moves away from the black hole forever, coming to rest at infinity. The collection of such particles provides an ambient medium with respect to which the motion of the jet clump can be described. That is, this ambient medium provides a dynamical way of characterizing the rest frame of the black hole.

Imagine two neighboring test particles at $`r`$ and $`r+\delta r`$ along the rotation axis of a central source. Far away from the source, the Newtonian gravitational acceleration of the reference particle at $`r`$ is $`GM/r^2`$ and for the other particle at $`r+\delta r`$ is $`GM/(r+\delta r)^2`$. The relative or tidal acceleration is then given by $`2GM\delta r/r^3`$, which by (10) is $`k\delta r`$ for $`r>>GM`$ and $`r>>a`$. Therefore, (7) and (10) reduce to the familiar Newtonian result once the relativistic tidal acceleration proportional to $`V^2<<1`$ is neglected.

The speed of a clump can be measured by observing the displacement of the clump in the flow relative to the ambient medium . Therefore, consider a clump along the axis of rotation moving rapidly past the reference particle belonging to the ambient medium. The equation of relative motion is given by (1) with

$$A_T=k(12V^2)Z,$$
(11)
where $`k`$ is determined via equations (9) and (10). Under our assumption we have that $`r^2>3a^2`$, and thus $`k<0`$. Also, we have assumed that $`\gamma =1`$ in (9) for the ambient medium surrounding the black hole; therefore, condition (8) is satisfied, the two special solutions of (11) are indeed attractors and the clump speed (predicted by our model) tends to $`0.7`$ over time . In  a similar assertion was also made for $`\gamma <1`$, which turns out to be erroneous in general; in any case, it is not relevant for a black hole source with $`aGM`$. The qualitative dependence of the solution of the generalized Jacobi equation upon $`\gamma 1`$ is depicted in Figure 2; in constructing this figure, we consider equation (11) simply as a nonlinear differential equation and ignore its possible physical limitations.

It is important to recognize two salient features of equation (11) for the force-free motion of an ultrarelativistic clump relative to the ambient medium with $`\gamma =1`$. First, the initial deceleration of the clump can be quite significant *within the range of validity* of the generalized Jacobi equation (11). For instance, Figure 3 depicts the deceleration of clumps with initial Lorentz factors of $`\mathrm{\Gamma }_0=20`$, 10 and 5, where we have introduced the Lorentz factor of the clump $`\mathrm{\Gamma }=(1V^2)^{1/2}`$, and $`\mathrm{\Gamma }_0`$ corresponds to the initial clump speed $`V_0`$ at $`r=r_0`$ and $`T=Z=0`$. Launched from $`r_0=10GM`$, the clumps decelerate to $`\mathrm{\Gamma }3`$ in about $`T=10GM5\times 10^5(M/M_{})\text{sec}`$. Although the distance of such a clump relative to the reference particle increases rapidly and soon the validity of the generalized Jacobi equation breaks down, the subsequent motion can be measured with the same model equation relative to the motion of another suitable reference particle in the ambient medium. This observation uses the second remarkable feature of equation(11): it does not explicitly depend on the choice of the reference particle in the ambient medium with $`\gamma =1`$. To illustrate this point, imagine a sequence of $`n`$ reference particles at $`r_i`$, $`i=1,2,\mathrm{},n`$, all with $`\gamma =1`$ along the rotation axis of the central source such that $`r_1<r_2<\mathrm{}<r_i<r_{i+1}<\mathrm{}<r_n`$. We can set up distinct Fermi frames along the worldline of each of these reference particles such that within each interval $`[r_i,r_{i+1}]`$ the generalized Jacobi equation is valid in the Fermi frame based on $`r_i`$. Once the clump reaches $`r_{i+1}`$, we can switch to the new Fermi frame based on $`r_{i+1}`$. So long as the initial speed of the clump relative to $`r_{i+1}`$ is above $`1/\sqrt{2}`$, the model predicts that the clump continues to decelerate toward the terminal speed. In this sense, the speed of the clump *relative to the ambient medium* approaches the terminal speed over time.

In the absence of any nongravitational forces, equations (1) and (11) express the geodesic equation for the center of mass of the clump in Fermi coordinates. Once the electromagnetic forces are included, the general equation of motion takes the form

$$\rho _m^{}(\frac{d^2x^\mu }{ds^2}+\mathrm{\Gamma }_{\nu \rho }^\mu \frac{dx^\nu }{ds}\frac{dx^\rho }{ds})=F_\sigma ^\mu j^\sigma ,$$
(12)
where $`\rho _m^{}`$ is the invariant mass density of the clump, $`F_{\mu \nu }`$ is the Faraday tensor of the exterior electromagnetic field in the ambient medium and $`j^\mu `$ is the charge current of the clump. Expressed in the Fermi coordinate system $`(T,𝐗)`$, equation (12) reduces to equation (1), where $`A_L`$ is due to the Lorentz force. In evaluating the Lorentz acceleration $`A_L`$ in the following section, we employ the physically reasonable approximation that the effects of spacetime curvature can be neglected for the sake of simplicity; thus, we treat the Fermi frame established along the reference geodesic as an inertial frame of reference in Section 3.

It is conceivable that the electromagnetic field configuration is such that $`F_\sigma ^\mu j^\sigma =0`$ in equation (12) and hence the motion of the clump is force-free (see, for example, ). This may be the case, for instance, beyond the launching point of the flow, a few gravitational radii from the central source. In the force-free case, the detailed analysis of  is applicable in the sense explained above and the clump speed tends to $`1/\sqrt{2}0.7`$, corresponding to a Lorentz factor of $`\sqrt{2}1.4`$. Our results may be of interest in connection with Galactic superluminal sources ; in fact, as already pointed out in , some of these superluminal jets have speeds that may be near $`1/\sqrt{2}0.7`$. In the rest of this paper, we investigate clump motion that is *not* force-free.

## 3 Lorentz Acceleration

Imagine a definite clump of plasma with constant mass $`m`$ moving with respect to the ambient medium with velocity $`𝐕=V\widehat{𝐧}`$, where $`\widehat{𝐧}`$ is a fixed direction in space that will be later identified with the rotation axis of the central gravitational source. The Lorentz force law for the motion of the clump can be written as

$$\frac{d}{dT}\left(\frac{mV}{\sqrt{1V^2}}\right)=(\delta 𝒱)(\rho 𝐄+𝐉\times 𝐁)\widehat{𝐧},$$
(13)
where $`\delta 𝒱`$, $`\rho `$ and $`𝐉`$ are, respectively, the volume, charge density and charge current of the clump with respect to the inertial frame $`(T,𝐗)`$. We assume that in the rest frame of the clump $`m=\rho _m^{}\delta 𝒱^{}`$, where $`\rho _m^{}`$ and $`\delta 𝒱^{}`$ are the proper density and volume of the plasma clump, respectively, and are assumed to be constants throughout the motion under consideration here. Primed quantities refer to the rest frame of the clump. It follows from Lorentz contraction that (13) may be written as

$$\frac{d}{dT}\left(\frac{V}{\sqrt{1V^2}}\right)=\frac{\sqrt{1V^2}}{\rho _m^{}}(\rho 𝐄+𝐉\times 𝐁)\widehat{𝐧}.$$
(14)
Moreover, we note that for one-dimensional motion

$$\frac{d}{dT}\left(\frac{V}{\sqrt{1V^2}}\right)=(1V^2)^{3/2}\frac{dV}{dT},$$
(15)
so that equation (14) reduces to

$$\frac{dV}{dT}=\frac{(1V^2)^2}{\rho _m^{}}(\rho 𝐄+𝐉\times 𝐁)\widehat{𝐧}.$$
(16)
We now turn to the electromagnetic aspects of (16) and assume that in the rest frame of the clump $`\rho ^{}=0`$ and that $`𝐉^{}=\sigma 𝐄^{}`$ in accordance with Ohm’s law. The electrical conductivity $`\sigma `$ is a tensor in general; however, we take $`\sigma `$ to be a scalar for the sake of simplicity. It follows from

$$\rho ^{}=\mathrm{\Gamma }(\rho 𝐕𝐉)$$
(17)
that $`\rho =𝐕𝐉=VJ_{}`$. This fact immediately implies that $`J_{}^{}=\mathrm{\Gamma }^1J_{}`$ using the transformation

$$J_{}^{}=\mathrm{\Gamma }(J_{}V\rho ).$$
(18)
Thus, we have

$$𝐉_{}^{}=𝐉_{},J_{}^{}=\mathrm{\Gamma }^1J_{}$$
(19)
and at the same time

$$𝐄_{}^{}=\mathrm{\Gamma }(𝐄+𝐕\times 𝐁)_{},E_{}^{}=E_{}.$$
(20)
It follows from (19), (20) and $`𝐉^{}=\sigma 𝐄^{}`$ that

$$𝐉=\sigma \mathrm{\Gamma }(𝐄+𝐕\times 𝐁).$$
(21)
Thus $`\rho =VJ_{}=\sigma \mathrm{\Gamma }VE_{}`$ and

$$(\rho 𝐄+𝐉\times 𝐁)\widehat{𝐧}=\sigma \mathrm{\Gamma }VE_{}^2+\sigma \mathrm{\Gamma }[(𝐄+𝐉\times 𝐁)\times 𝐁]_{}.$$
(22)
Using the identity

$$[(𝐕\times 𝐁)\times 𝐁]_{}=V(B^2B_{}^2)=VB_{}^2,$$
(23)
we find that

$$(\rho 𝐄+𝐉\times 𝐁)\widehat{𝐧}=\mathrm{\Gamma }\sigma B^2(w\frac{B_{}^2E_{}^2}{B^2}V),$$
(24)
where

$$w=\frac{(𝐄\times 𝐁)\widehat{𝐧}}{B^2}$$
(25)
and

$$𝐰_D=\frac{𝐄\times 𝐁}{B^2}$$
(26)
is the electrical drift velocity . To proceed further, we need a model for the electromagnetic field configuration in the medium surrounding the black hole. This is considered in the next section.

## 4 A Simple Model

Let us now consider the magnetically dominated environment in which the clump is moving. In some models for such an environment  and in cylindrical coordinates $`(\rho ,\varphi ,z)`$ with $`\widehat{𝐧}=\widehat{𝐳}`$, integrating Maxwell’s equation $`𝐁=0`$ implies that $`B_z\rho ^2`$. If angular momentum is approximately conserved in the electromagnetic field, then $`B_\varphi \rho ^1`$. Thus, at large distances from the source, jets and collimated winds are often considered to have predominantly toroidal magnetic fields. Assuming that the plasma is efficient at maintaining a very small proper electric field through its conductive properties, this implies that $`𝐄`$ is approximately radial in cylindrical coordinates . It is clear from this brief description that in any reasonably realistic scenario the electromagnetic field configuration would be rather complicated. However, we seek a simple model situation that would render the resulting equation of motion amenable to mathematical analysis. Therefore, we assume that the average electric field is primarily radial, $`𝐄=E\widehat{𝝆}`$, and the average magnetic field is primarily azimuthal, $`𝐁=B\widehat{\mathit{\varphi }}`$, such that $`E<B`$. In view of the symmetry of the central source about its equatorial plane, in this paper we concentrate—for the sake of simplicity—only on the clump moving along the positive $`z`$-axis. It follows from (25) that the electrical drift speed is given by $`w=E/B<1`$. Moreover, (24) takes the simple form $`\mathrm{\Gamma }(\sigma B^2)(wV)`$.

It turns out that in the rest frame of the clump $`𝐄^{}=E^{}\widehat{𝝆}`$ and $`𝐁^{}=B^{}\widehat{\mathit{\varphi }}`$, i.e. the fields have the same configuration as in the global inertial frame $`(T,𝐗)`$. Thus following the analysis in §2.10 of , the conductivity $`\sigma `$ should be identified with $`\sigma _{}`$ as the current $`𝐉^{}=\sigma _{}𝐄^{}`$ will be crossing the magnetic lines of force in the local rest frame of the clump. On the other hand, for motion along the magnetic lines of force, the electrical conductivity $`\sigma _{}`$ would have its usual value

$$\sigma _{}=\frac{n_ee^2\tau _c}{\mu _e},$$
(27)
where $`e`$, $`n_e`$ and $`\mu _e`$ are, respectively, the charge, number density and mass of the electron in the electron-proton plasma in the rest frame of the clump. Here the relaxation time $`\tau _c=\nu _c^1`$ is the average time interval between electron collisions and hence $`\nu _c`$ is the average frequency of electron collisions per second. It follows from the discussion in §2.10 of  that in case of a tenuous plasma in a strong magnetic field such that $`\nu _c^2\mathrm{\Omega }_e^2`$, we have

$$\frac{\sigma _{}}{\sigma _{}}=\left(\frac{\nu _c}{\mathrm{\Omega }_e}\right)^2,$$
(28)
where $`\mathrm{\Omega }_e=eB^{}/\mu _e`$ is the cyclotron frequency for a free electron. Thus

$$\frac{\sigma _{}B^2}{\rho _m^{}}=\left(\frac{\mu _e}{\mu _p}\right)\left(\frac{B}{B^{}}\right)^2\nu _c,$$
(29)
where $`\mu _p`$ is the proton mass and $`\rho _m^{}=n_p\mu _p+n_e\mu _e=n_e(\mu _p+\mu _e)n_e\mu _p`$, since $`\rho ^{}=0`$ implies that the density of protons in the clump is the same as that of electrons, $`n_p=n_e`$, and $`\mu _p/\mu _e18361`$. Moreover,

$$𝐁_{}^{}=\mathrm{\Gamma }(𝐁𝐕\times 𝐄)_{},$$
(30)
which implies that

$$\frac{B^{}}{B}=\mathrm{\Gamma }(1wV).$$
(31)
Combining these results with (16) reveals that the Lorentz acceleration is given by

$$A_L=\frac{dV}{dT}=\alpha \frac{(1V^2)^{5/2}}{(1wV)^2}(wV),$$
(32)
where

$$\alpha =\left(\frac{\mu _e}{\mu _p}\right)\nu _c.$$
(33)
It is interesting to note that for copper at room temperature $`\nu _c10^{14}\text{sec}^1`$.

It follows from equations (32) and (33) that the Lorentz acceleration $`A_L`$ vanishes if the clump speed is equal to the drift speed. Moreover, $`A_L`$ is proportional to the collision frequency $`\nu _c`$. This is due to the fact that only through collisions can an appropriate interior current be established that could lead to the acceleration of the clump along the flow direction.

It remains to provide an estimate of $`\nu _c`$ for the case of the electron-proton plasma under consideration here. This is a rather complicated problem and we therefore limit our investigation to the simpler case of collisions among the electrons while neglecting the motion of the protons. Thus we assume that $`\tau _c`$ can be estimated using the “self-collision time” of electrons given by equation (5-26) of

$$\tau _c^{}=\frac{0.266𝒯^{3/2}}{n_e\mathrm{ln}\mathrm{\Lambda }}\text{sec},$$
(34)
where $`𝒯`$ is the absolute temperature and $`n_e`$ is the number of electrons in the clump per $`\text{cm}^3`$. Here $`\mathrm{\Lambda }`$ is essentially the ratio of two lengths: the Debye shielding distance over the characteristic impact parameter of electron collisions such that the deflection angle in the orbital plane is equal to $`\pi /2`$. One can use Table 5.1 of  to find the values of $`\mathrm{ln}\mathrm{\Lambda }`$ in terms of $`𝒯`$ and $`n_e`$. According to this table, for a fixed $`𝒯`$, $`\mathrm{ln}\mathrm{\Lambda }`$ decreases very slowly but monotonically as $`n_e`$ increases from $`n_e=1\text{cm}^3`$ to $`n_e=10^{24}\text{cm}^3`$; for example, for $`𝒯=10^8\text{K}`$, $`\mathrm{ln}\mathrm{\Lambda }`$ monotonically decreases from $`34.3`$ for $`n_e=1\text{cm}^3`$ to $`6.69`$ for $`n_e=10^{24}\text{cm}^3`$. Thus for $`𝒯=10^8\text{K}`$ and $`n_e=10^6\text{cm}^3`$, we find $`\mathrm{ln}\mathrm{\Lambda }=27.4`$ and hence from (34) with $`\nu _c(\tau _c^{})^1`$ we obtain $`\alpha 0.56\times 10^7\text{sec}^1`$. For fixed $`𝒯`$, $`\alpha `$ increases almost linearly with $`n_e`$ such that for $`n_e=10^{15}\text{cm}^3`$ in the case under discussion, we have $`\alpha 35\text{sec}^1`$ . These considerations illustrate the fact that $`\alpha `$ could have a considerable range of values depending upon the electron temperature $`𝒯`$ and density $`n_e`$.

## 5 Equation of Motion

The equation of motion of the clump is given by

$$\frac{d^{\mathrm{\hspace{0.17em}2}}Z}{dT^2}+k(T)(12V^2)Z=\alpha \frac{(1V^2)^{5/2}}{(1wV)^2}(wV),$$
(35)
where $`V=dZ/dT`$, $`\alpha `$ is a characteristic of the plasma clump, $`w`$ is a characteristic of the electromagnetic field of the ambient medium and $`k`$ is essentially the curvature of the central source given by

$$k(T)=\frac{2GMz(z^23a^2)}{(z^2+a^2)^3},$$
(36)
while $`z(T)`$ is determined using the geodesic equation for the background

$$\left(\frac{dz}{dT}\right)^2=\gamma ^21+\frac{2GMz}{z^2+a^2}.$$
(37)
The equation of motion (35) has been derived for the clump moving along the positive $`Z`$-axis. Recall that $`z`$ in equations (36) and (37) is a substitute for the radial coordinate $`r`$. Therefore, let us note that equations (35)–(37) are invariant under the transformations $`ZZ`$, $`ww`$, $`VV`$ and $`zz`$. Hence our results hold for the clump moving along the negative $`Z`$-direction as well.

It is useful to put equations (35)–(37) in dimensionless form. To this end, we assume that all lengths are measured in units of $`GM`$, which is one-half of the gravitational radius of the central source. Defining dimensionless quantities $`\widehat{k}`$, $`\widehat{\alpha }`$ and $`\widehat{a}`$ as

$$\widehat{k}=(GM)^2k,\widehat{\alpha }=GM\alpha ,\widehat{a}=\frac{a}{GM}$$
(38)
we recover equations (35)–(37) in dimensionless form once we let $`k\widehat{k}`$, $`\alpha \widehat{\alpha }`$, $`a\widehat{a}`$ and $`GM1`$.

The equations of motion in dimensionless form are equivalent to the dynamical system

$`{\displaystyle \frac{dz}{dT}}`$ $`=`$ $`\left(\gamma ^21+{\displaystyle \frac{2z}{z^2+\widehat{a}^2}}\right)^{1/2},`$
$`{\displaystyle \frac{dZ}{dT}}`$ $`=`$ $`V,`$
$`{\displaystyle \frac{dV}{dT}}`$ $`=`$ $`\widehat{\alpha }{\displaystyle \frac{(1V^2)^{5/2}}{(1wV)^2}}(wV)+{\displaystyle \frac{2z(z^23\widehat{a}^2)}{(z^2+\widehat{a}^2)^3}}(12V^2)Z,`$ (39)
where we assume that $`\gamma 1`$.

We have proved the following result: If $`\widehat{\alpha }>0`$, $`0w<1`$, $`z(0)>\sqrt{3}\widehat{a}`$, and $`0<V(0)<1`$, then $`lim_T\mathrm{}V(T)=w`$. That is, as time increases the clump speed approaches the drift speed of the ambient medium.

Equation (35) and the system (5) show how the tidal equation can in principle be generalized to include nongravitational forces. On the other hand, it is necessary to emphasize the qualitative significance of equation (35), since its gravitational part is valid only within a distance $``$ of the fiducial particle and its electromagnetic part is based on a rather simple model. If the tidal part is completely ignored, then the Lorentz force would lead to a monotonic deceleration of an initially ultrarelativistic clump toward $`w`$. Numerical experiments based on the system (5) demonstrate that the tidal part initially dominates, causing a very rapid drop in $`\mathrm{\Gamma }`$ toward $`\sqrt{2}`$, and this occurs mostly within the domain of validity of the tidal term, but then the electromagnetic term takes over and $`\mathrm{\Gamma }`$ slowly tends toward $`(1w^2)^{1/2}`$. Figure 4 depicts a profile for the graph of $`\mathrm{\Gamma }`$ versus time, where the parameter values and initial conditions are chosen to be typical for a microquasar with $`M=10M_{}`$. According to our model (5) with electrical drift speed corresponding to $`\mathrm{\Gamma }=2`$ and with Lorentz acceleration coefficient $`\widehat{\alpha }=10^7`$, an initial Lorentz factor of $`\mathrm{\Gamma }=5`$ decreases to $`\mathrm{\Gamma }1.5`$ in about $`0.25\text{sec}`$ and then very slowly increases to within 1% of $`\mathrm{\Gamma }=2`$ after 48 hours. More generally, for $`\widehat{\alpha }<10^2`$ the timescale for the initial rapid decrease in $`\mathrm{\Gamma }`$ is about $`0.1\text{sec}`$. The corresponding total relaxation time to $`\mathrm{\Gamma }=2`$ is approximately $`(50\widehat{\alpha })^1\text{sec}`$, so that for $`\widehat{\alpha }`$ ranging from $`10^3`$ to $`10^{11}`$, it ranges from $`20\text{seconds}`$ to $`64\text{years}`$. For a quasar with $`M=10^8M_{}`$, the timescale for the initial rapid decrease in $`\mathrm{\Gamma }`$ would be about two weeks.

To see that our result holds in general, let us note that $`dz/dT>0`$ for $`T>0`$. Thus, $`z`$ is an increasing function of $`T`$, which is defined for $`0T<\mathrm{}`$. By viewing $`z`$ as a new “time” in the dynamical system (5), it suffices to consider the system

$`{\displaystyle \frac{dZ}{dz}}`$ $`=`$ $`{\displaystyle \frac{1}{R(z)}}V,`$
$`{\displaystyle \frac{dV}{dz}}`$ $`=`$ $`{\displaystyle \frac{1}{R(z)}}\left(\widehat{\alpha }{\displaystyle \frac{(1V^2)^{5/2}}{(1wV)^2}}(wV)+{\displaystyle \frac{2z(z^23\widehat{a}^2)}{(z^2+\widehat{a}^2)^3}}(12V^2)Z\right),`$ (40)
where

$$R(z):=\left(\gamma ^21+\frac{2z}{z^2+\widehat{a}^2}\right)^{1/2}.$$

In this system $`dZ/dz>0`$. Hence, as before, it suffices to consider the scalar differential equation

$$\frac{dV}{dZ}=\widehat{\alpha }\frac{(1V^2)^{5/2}}{(1wV)^2V}(wV)+2Z\frac{S(Z)(S(Z)^23\widehat{a}^2)(12V^2)}{(S(Z)^2+\widehat{a}^2)^3V},$$
(41)
where $`S`$ is the function such that $`S(Z(z))=z`$.

Direction fields of $`dV/dZ`$ in the $`(Z,V)`$-plane are depicted in Figure 5 for the two cases $`w>1/\sqrt{2}`$ and $`w<1/\sqrt{2}`$. Graphs of typical solutions are also drawn. Since $`V`$ stays bounded and the curvature $`k(T)`$ approaches zero very fast, it can be shown that the second term on the right-hand side of equation (41) rapidly approaches zero as $`T\mathrm{}`$. Also, with this term set to zero, the dynamical system (41) has an asymptotically stable steady state at $`V=w`$. For large $`Z`$, the first term is dominant; therefore, the limiting value of the solution is $`w`$.

## 6 Discussion

This paper is devoted to a simple model of the evolution of the speed of a clump in a relativistic flow once it emerges from the environment immediately surrounding the central source. Combined with a realistic treatment of plasma effects, our theoretical approach should be of interest in the study of astrophysical jets. If the motion of the clump is force-free, then it follows from the repeated application of the generalized Jacobi equation that the clump speed tends to $`V_T=1/\sqrt{2}`$ corresponding to $`\mathrm{\Gamma }=\sqrt{2}`$. On the other hand, if the force-free condition is not satisfied, the clump speed tends to $`V_L=w`$, where $`w`$ is the electrical drift speed. In fact, if the tidal component initially dominates, as would be the case for most microquasars, then starting with an initial speed almost equal to unity, the tidal force is responsible for the initial rapid decrease of the velocity toward $`1/\sqrt{2}`$, but then the Lorentz term takes over and the velocity slowly approaches $`w`$. On the other hand, if the Lorentz term is dominant, as would be the case for most quasars, there is in effect a relatively slow decrease of the velocity toward $`w`$.


format: text

text:
# Radiation environment along the INTEGRALBased on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), Czech Republic and Poland, and with the participation of Russia and the USA. orbit measured with the IREM monitor

## 1 Introduction

The INTEGRAL Radiation Environment Monitor (IREM) is a space dedicated detector assembly for both dosimetry and on board electron and proton spectroscopy. It is an adapted version of the Standard Radiation Environment Monitor (SREM) developed in partnership between European Space Agency (ESA), PSI and Contraves Space AG (Zürich) 1 (Contraves, 1996). Ten identical instruments were manufactured and calibrated for the ESA space program. Three monitors are by now in space (on board STRV, on board PROBA and on board INTEGRAL) and the rest have already been assigned to the forthcoming scientific missions. As it is well illustrated by XMM-Newton, CHANDRA or INTEGRAL itself, an autonomous radiation monitoring is of great importance for the spacecraft operations as well as for the lifetime and health of its instruments and devices. Having a reliable radiation monitor on board can help optimizing the mission observing time, react quickly to elevated radiation levels (by providing indications as to when safety measures must be applied for the most sensitive devices) or support in tracing the spacecraft anomalies.

The INTEGRAL 9 (Winkler, 2003) mission’s prime goals are studies of intense gamma radiation sources and explorations of rare and powerful events. Four very sensitive payload instruments allow for observations in the range 3 keV to 10 MeV. In addition, an optical monitor (500-850 nm) allows observation in the V band. In order to maximize uninterrupted observing time and protect vulnerable equipment from hazardous radiation the satellite’s orbit extends from 10 000 to 153 000 km. It allows spending almost 90% of its 72 hours long revolution outside of the Earth’s radiation belts (altitude above 40 000 km). The IREM is responsible for the radiation monitoring on board and its key function on INTEGRAL is a continuous checking of the radiation environment to alert the spacecraft when high radiation levels are met. The payload instruments rely on such information and react accordingly entering, if necessary, a special safe mode in which they are protected from possible radiation damages.

The IREM was switched on only 10 hours after the launch (2002 October the 17th) and after a short commissioning phase it began its routine operation. During the INTEGRAL mission it can achieve its two objectives: being a vital part of the spacecraft radiation protection system and functioning as an autonomous radiation monitoring device. The mission highly elliptical orbit allows IREM to probe both the dynamic outer electron belt and the interplanetary environment where cosmic rays, solar proton and electron events (as well as other phenomena like energetic Jovian electrons) are encountered. With the planned mission lifetime of up to 5 years the IREM will be able to cover the whole declining phase of the current Solar cycle.

We will briefly characterize below the main features of the monitor, including its calibration and response modeling 2 (Hajdas et al., 2002) and present its first measurements of the external radiation environment along the INTEGRAL orbit.

## 2 Instrument characteristics

The IREM (Figure 1) utilizes three standard Silicon Surface Barrier Detectors; 500 $`\mu `$m thickness, active area 25 (2) and 50 (1) mm2 \- see IREM Users Manual111e.g. http://pif.web.psi.ch/irem/IREM\_USER\_MANUAL\_v3.pdf for more details. They are embedded in a bi-metallic shielding of Tantalum (inner) and Aluminum (outer) with 8 g/cm2 of total thickness. For enhanced resolution in energy and directionality of the incident particles, two of the detectors are arranged in a telescope. All pre-amplified detector pulses are scrutinized by a set of fifteen comparators - ten for single events, four for coincidences and one heavy ion channel. Their levels are optimized to get the most accurate information on the spectral shape of the detected particles. The low energy detection thresholds are: $`E_{\mathrm{thr}}^\mathrm{p}`$ 10 MeV for protons and $`E_{\mathrm{thr}}^\mathrm{e}`$ 0.5 MeV for electrons. The heavy ion channel has an energy threshold of $``$ 150 MeV/nucleon. The particles come through the conical front collimators of $`\pm 20\mathrm{°}`$ opening.
High energy particles ($`E_\mathrm{p}>`$100 MeV) can enter the detector from any direction. Some of them, however, are stopped in the satellite bulk mass before they hit the monitor. Therefore, the full response matrix must take the satellite into account 2 (Hajdas et al., 2002).

## 3 Calibrations and modeling

IREM calibrations were done using the Proton Irradiation Facility (PIF) 4 (Hajdas et al. 1996) as well as gamma and electron radioactive sources in PSI. It was important to use the same particles and spectra as anticipated during the mission in space. The full data set comprising several initial particle energies ($`E_0`$: 8-300 MeV) and incoming angles ($`\theta `$: 0-180$`\mathrm{°}`$) provided a reference for the response matrix. Experimental results were compared with fine-tuned computer calculations performed with the help of the GEANT code from CERN using a precise computer model of the monitor. In the next step, IREM responses for both protons and electrons were generated for the whole energy range anticipated in space.

The proton response was calculated for thirty energy bins equally spanned on logarithmic scale between 8 and 800 MeV. The response presented for selected scalers in Figure 2 (lower panel) shows results integrated over the full $`4\pi `$ angle of incoming particles. The right hand side of the sensitivity curve strongly depends on the extra shielding provided by the satellite. For electrons, the number of bins was only equal to 15 and covered the energy range from 0.3 to 15 MeV. The bins had equal widths on logarithmic scale. Response calculations (also integrated over the full $`4\pi `$ angle) are presented for selected scalers in Figure 2 (upper panel). The shielding of the monitor stops not only the bremsstrahlung but also electrons coming from outside of the entrance collimators. It implies that their response function is only very weakly affected by the satellite.

## 4 IREM performance and alerts

The IREM was the first instrument on board to be switched on after the launch. Its basic task of warning the spacecraft payload instruments when high levels of radiation are met is realized by sending periodic broadcast packets. Every eight seconds the IREM passes to the spacecraft a fifteen words long Transfer Data Block (TDB). Its first five words contain information about the current radiation environment as well as about the IREM status. The rest consists of either the scientific or housekeeping data for further download to the ground. Three TDB words containing radiation data come from pre-selected, dead time corrected IREM scalers. One of them monitors proton flux, the other one is sensitive to electrons, while the third one informs about a deep dose deposition. Each payload instrument has its own response method depending on individual radiation hardness - see Table 1.

In addition to radiation warning signals from IREM the mission ground station provides a Radiation Belt alert flag (currently set for descending direction - belt entrance - at 60 000 km and ascending - belt exit - at 40 000 km) and all payload devices have an independent safeguard too.

## 5 Radiation environment

The orbit of INTEGRAL was selected to maximize its uninterrupted scientific observing time and telemetry flow, reducing at the same time its encounter with the Earth’s radiation belts to minimum. High energy particles in the belts cause not only an enhanced instrument background but also induce radiation damages and malfunctioning of the spacecraft devices. It may strongly diminish an effective lifetime and data quality of the space observatory. The satellite revolution of 72 hours duration has a perigee of about 10 000 km. The trajectory crosses the whole outer electron belt and just barely touches the inner belt with increased fluxes of protons.

As one can see in Figure 3, the increase of the count rate inside of the belts may reach 5 orders of magnitude and it is caused by high energy electrons. With an apogee of 153 000 km the mission spends up to 90% of the orbital time in regions dominated only by a fairly constant cosmic rays background (mostly high energy protons with a flux of $``$ 2 cm-2s-1). The IREM used cosmic rays as a cross-check of its detector calibration. The mission observing program is conducted outside the belts. It may, however be interrupted by infrequent solar events.

## 6 Radiation belts

The spacecraft passage through the belts is characterized by highly variable radiation environment 5 (Bühler & Desorgher, 2002) that depends upon geomagnetic coordinates of each particular trajectory. A typical belt crossing takes between seven and ten hours and electron fluxes may reach up to 107 cm-2s-1 (for energies higher than 0.5 MeV). The maximum flux shows variations from orbit to orbit by a factor of about three. Due to the high perigee value, the IREM observes only a slight increase of the proton flux from the outer edge of the proton belt. Its intensity reaches just about $`10`$ cm-2s-1 for energies above 20 MeV. It takes about 1.5 hours to fly through the proton populated area.

The standard coordinate system used to display or analyze radiation belt particle fluxes is making use of the 2 following variables: the magnetic field strength $`B`$ and the McIlwain $`L`$-shell parameter $`L`$ 6 (McIlwain, 1966). For a given satellite position, parameters $`B`$ and $`L`$ are computed using the International Geographic Reference Field, IGRF plus an external field model (representing the solar wind influenced parts of the Earth’s magnetic field). The IREM electron spectra are usually approximated by an exponential function. The equation below is used to describe the differential electron flux $`f`$ (cm-2s-1MeV-1) as a function of energy $`E`$

$$f(E)=Ne^{\gamma (E1.0\mathrm{MeV})}.$$
(1)
Parameters $`N`$ and $`\gamma `$ are determined from the measured count rates by a fitting procedure. They depend on the IREM detection system sensitivity and allow to parameterise the electron spectra between 0.5 and 5 MeV. In Figure 4, the average electron flux normalization parameter $`N`$ and the spectrum hardness parameter $`\gamma `$ are shown as function of $`L`$-values, for INTEGRAL passages through the outer radiation belt. The data illustrates a period from January 1 to May 1 2003. The spectra are hardest at low $`L`$-values i.e. closer to the Earth and soften with increasing values of $`L`$ while the maximum fluxes are found for $`L`$ $``$ 4-5.

Monitoring of the electron belts boundaries is one of the most important functions of the IREM. In addition to protecting sensitive devices it allows for gaining extra observing slots (up to 3% of observing time per orbit) as the scientific program of the mission is conducted only outside of the belts. Time evolution of the belt entrance and exit limits for the first 81 spacecraft orbits is presented in Figure 5. Large deviations around the mean value are attributed to both random (like solar flares) and periodic (like its rotation) mechanisms of the solar activity as well as to the sun-earth geomagnetic bond with its seasonal variations.

The radiation environment in space, as measured using on board monitors like the IREM, can be compared with existing models of the belts. The electron AE-8 7 (Sawyer & Vette, 1976) and proton AP-8 8 (Vette, 1999) NASA belt models are quasi-standards, conventionally used to asses the radiation environment on spacecrafts. They both are static representations of average fluxes of particles trapped in the Earth’s radiation belts.

For the proton belt the static approximation is usually qualitatively good enough but this may not be the case for the electron belts as they are highly dynamic 5 (Bühler & Desorgher, 2002). This is demonstrated in Figure 6 where the IREM measured count rates are compared with predictions computed with the AP8/AE8 models. The upper figure shows count rates in the TC3 electron counter (low energy threshold of $``$ 500 keV) for the perigee passage of May 4 2003. Although the measured particle rates are much higher than the predicted ones, both the peak positions and distributions are similar. The lower figure shows an analogous passage (31 May 2003) after a magnetic storm that was initiated by a solar proton event. In this case, neither intensities or peaks nor distributions of measured and predicted environment agree.

## 7 Solar events

The time spent outside of the magnetosphere is dominated by cosmic rays and may be occasionally interrupted by particles ejected from the sun during solar events. It is represented by a flat region between the belts peaks as seen in Figure 3. High energy particle fluxes from coronal mass ejections may have duration from hours to days and can heavily disturb the observation program. Several events of different amplitude and duration have already occurred during the INTEGRAL mission reflecting the fact that the Sun is still quite active (maximum part of the Solar cycle). Their long term impact on the radiation belts was already illustrated in Figure 6 while a light curve of one such solar event itself is shown in Figure 7 in which one can see an increase of the count rates by more than 100% in scalers TC3 and S14.

In the upper panel of Figure 7, one could see the contribution of low energy protons ($`E_\mathrm{p}`$$`>`$10 MeV) to the TC3 count rate (mainly sensitive to electrons $`E_\mathrm{e}`$$`>`$500 keV) while the middle panel present the higher energy ($`E_\mathrm{p}`$$`>`$20 MeV) proton count rate from S14. The peak on the left side of both graphs is due to radiation belt particles - electrons in case of TC3 and protons for S14 while the long structure in the middle is made by solar protons. The bottom panel shows the counters ratio in which one sees the evolution of the proton hardness. The peak in the hardness occurred about 5 hours before arrival of the maximum flux.

## 8 Other instruments

The IREM calibration of count rates was confirmed by other payload instruments during commissioning phase of the mission. For this purpose one used the active shielding BGO detectors of the high energy instruments as well as the CCD chips of the OMC. The verification was performed with highly penetrating cosmic rays and radiation belt particles. Further cross checks are routinely performed using for example the JEM-X and ISGRI data while entering the radiation belt, when the devices are still in their active modes. Such verification revealed very stable performance of all IREM detectors.

## 9 Summary

The IREM, flying on board INTEGRAL, belongs to the ESA Standard SREM monitors that are optimized for detection of protons and electrons and for alerting the spacecraft during high radiation levels. It is specially suited to the INTEGRAL payload radiation protection scheme. The IREM permanently monitors an external radiation environment of the satellite and periodically sends broadcast packets with current levels of particle fluxes. Each payload instrument reacts to the IREM message individually.
Most of the spacecraft orbit is characterized by a quasi constant cosmic rays background that may be occasionally disturbed by sporadic CME solar events. When approaching a perigee, the spacecraft passes through the outer electron belt. IREM measurements reveal very dynamic belt environment that shows only a qualitative agreement with the present NASA AP8/AE8 models. Real-time IREM radiation maps allow not only for scientific program optimizing and instrument protection. The data are also used for the space weather global programs and are promptly (i.e. within 2 hours) available for further analysis.


text:
# 1 Introduction

## 1 Introduction

Multipole expansion of an incident radiation field in terms of vector spherical wavefunctions (VSWFs) is required for theoretical scattering methods such as the T-matrix method  and generalised Lorenz-Mie theory (GLMT) . Since the VSWFs are a complete orthogonal basis for solutions of the vector Helmholtz equation,

$$^2𝐄+k^2𝐄=0,$$
(1)
it is theoretically straightforward to find the multipole expansion for any monochromatic radiation field using the orthogonal eigenfunction transform , also known as the generalised Fourier transform. For example, the usual formula for multipole expansion of a plane wave can be derived in this way.

While the case of plane wave illumination is useful for a wide range of scattering problems, many applications of scattering involve laser beams. In particular, laser trapping  requires strongly focussed beams. As scattering calculations allow the optical forces and torques within the optical trap to be determined , and the T-matrix method is particularly useful for the repeated calculations typically required, multipole expansion of laser beams, including strongly focussed beams, is extremely useful.

Unfortunately, laser beams present some serious theoretical difficulties. These problems are not a result of any deficiency in the basic process of multipole expansion of radiation fields, but are due to the fact that standard representations of laser beams are not radiation fields—that is, their standard mathematical forms are not solutions of the vector Helmholtz equation, but are solutions of the paraxial scalar wave equation (higher order corrections can be used to improve the accuracy as reality becomes less paraxial ). Such pseudo-fields are only approximate solutions of the vector wave equation. The deviation from correctness increases as the beam is more strongly focussed .

While, in principle, any radiation field can be expanded in terms of multipoles, multipole expansions do not exist for approximate pseudo-fields that do not satisfy the vector wave equation. Since the standard paraxial representations of laser beams are so widely used, multipole expansions that correspond to the standard beams in a meaningful way are highly desirable. For this to be possible, some method must be used to “approximate” the standard laser beam with a real radiation field. We can note that there exists a significant and useful body of work on multipole expansion of beams . However, while satisfactorily efficient and accurate methods exist for weakly focussed beams, when the deviations from paraxiality are small, strongly focussed beams, for example, as required for laser trapping, remain problematic.

The traditional method used in scattering calculations was to assume that the actual incident field was equal to the paraxial field on the surface of the scatterer. This has the unfortunate drawback that the multipole expansion coefficients depend on the size, shape, location, and orientation of the scattering particle. As a result, an artificial dependence on particle position and size is introduced into scattering calculations and calculations of optical forces, since the differing multipole expansions correspond to different beams. This is obviously undesirable. While matching the fields on the surface of the scattering particle is not an adequate solution to the problem, the basic concept—matching the fields on a surface—is sound. However, a scatterer-independent surface must be used. Two natural choices present themselves: the focal plane (for beams with a well-defined focal plane), and a spherical surface in the far field.

Integral methods can be used over these surfaces (direct application of the orthogonal eigenfunction transform). Here, however, we use a point-matching method since, firstly, the method tolerates a much coarser grid of points, thereby increasing computational efficiency, and secondly, increased robustness can be gained by using extra points to give an overdetermined system of equations which can then be solved in a least-squares sense. We use both focal plane matching and far field matching, and evaluate their respective merits.

## 2 Point-matching method

All monochromatic radiation fields satisfy the vector Helmholtz equation (1) in source-free regions. The vector spherical wavefunctions (sometimes called vector spherical harmonics) are a complete set of orthogonal solutions to this equation. The VSWFs are

$`𝐌_{nm}^{(1,2)}(k\mathrm{r})`$ $`=`$ $`N_nh_n^{(1,2)}(kr)𝐂_{nm}(\theta ,\varphi )`$ (2)
$`𝐍_{nm}^{(1,2)}(k\mathrm{r})`$ $`=`$ $`{\displaystyle \frac{h_n^{(1,2)}(kr)}{krN_n}}𝐏_{nm}(\theta ,\varphi )+`$ (3)
$`N_n\left(h_{n1}^{(1,2)}(kr){\displaystyle \frac{nh_n^{(1,2)}(kr)}{kr}}\right)𝐁_{nm}(\theta ,\varphi )`$
where $`h_n^{(1,2)}(kr)`$ are spherical Hankel functions of the first and second kind, $`N_n=1/\sqrt{n(n+1)}`$ are normalisation constants, and $`𝐁_{nm}(\theta ,\varphi )`$, $`𝐂_{nm}(\theta ,\varphi )`$, and $`𝐏_{nm}(\theta ,\varphi )`$ are the vector spherical harmonics:

$`𝐁_{nm}(\theta ,\varphi )`$ $`=`$ $`𝐫Y_n^m(\theta ,\varphi )=\times 𝐂_{nm}(\theta ,\varphi )`$ (4)
$`=`$ $`\widehat{\theta }{\displaystyle \frac{}{\theta }}Y_n^m(\theta ,\varphi )+\widehat{\varphi }{\displaystyle \frac{\mathrm{i}m}{\mathrm{sin}\theta }}Y_n^m(\theta ,\varphi ),`$
$`𝐂_{nm}(\theta ,\varphi )`$ $`=`$ $`\times \left(𝐫Y_n^m(\theta ,\varphi )\right)`$ (5)
$`=`$ $`\widehat{\theta }{\displaystyle \frac{\mathrm{i}m}{\mathrm{sin}\theta }}Y_n^m(\theta ,\varphi )\widehat{\varphi }{\displaystyle \frac{}{\theta }}Y_n^m(\theta ,\varphi ),`$
$`𝐏_{nm}(\theta ,\varphi )`$ $`=`$ $`\widehat{𝐫}Y_n^m(\theta ,\varphi ),`$ (6)
where $`Y_n^m(\theta ,\varphi )`$ are normalised scalar spherical harmonics. The usual polar spherical coordinates are used, where $`\theta `$ is the co-latitude measured from the $`+z`$ axis, and $`\varphi `$ is the azimuth, measured from the $`+x`$ axis towards the $`+y`$ axis. The $`h_{n1}^{(1,2)}(kr)nh_n^{(1,2)}(kr)/kr`$ term in (3) results from the identity $`(\mathrm{d}/\mathrm{d}kr)krh_n^{(1,2)}(kr)=krh_{n1}^{(1,2)}(kr)nh_n^{(1,2)}(kr)`$.

$`𝐌_{nm}^{(1)}`$ and $`𝐍_{nm}^{(1)}`$ are outward-propagating TE and TM multipole fields, while $`𝐌_{nm}^{(2)}`$ and $`𝐍_{nm}^{(2)}`$ are the corresponding inward-propagating multipole fields. Since these wavefunctions are purely incoming and purely outgoing, each has a singularity at the origin. Since the incident field has equal incoming and outgoing parts, and must be singularity free in any source-free region, it is useful to define the singularity-free regular vector spherical wavefunctions:

$`\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r})`$ $`=`$ $`\frac{1}{2}[𝐌_{nm}^{(1)}(k\mathrm{r})+𝐌_{nm}^{(2)}(k\mathrm{r})],`$ (7)
$`\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r})`$ $`=`$ $`\frac{1}{2}[𝐍_{nm}^{(1)}(k\mathrm{r})+𝐌_{nm}^{(2)}(k\mathrm{r})].`$ (8)
Since the spherical Bessel functions $`j_n(kr)=\frac{1}{2}(h_n^{(1)}(kr)+h_n^{(2)}(kr))`$, the regular VSWFs are identical to the incoming and outgoing VSWFs except for the replacement of the spherical Hankel functions by spherical Bessel functions.

Since the VSWFs are a complete set of solutions to the vector Helmholtz equation (1), every solution can be written as a linear combination of VSWFs. Thus, the incident field can be written in terms of *expansion coefficients* (also called *beam shape coefficients*) $`a_{nm}^{(3)}`$ and $`b_{nm}^{(3)}`$ as

$$𝐄_{\mathrm{inc}}(\mathrm{r})=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=n}{\overset{n}{}}a_{nm}^{(3)}\mathrm{𝐑𝐠𝐌}_{nm}(k\mathrm{r})+b_{nm}^{(3)}\mathrm{𝐑𝐠𝐍}_{nm}(k\mathrm{r}).$$
(9)
Alternatively, since the outgoing field is purely a consequence of the incoming field (with or without a scatterer present), all the necessary information can be conveyed by writing the incident field expansion in terms of incoming wavefunctions as

$$𝐄_{\mathrm{inc}}(\mathrm{r})=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=n}{\overset{n}{}}a_{nm}^{(2)}𝐌_{nm}^{(2)}(k\mathrm{r})+b_{nm}^{(2)}𝐍_{nm}^{(2)}(k\mathrm{r}).$$
(10)
The two sets of expansion coefficients are related, since $`a_{nm}^{(3)}=2a_{nm}^{(2)}`$ and $`b_{nm}^{(3)}=2b_{nm}^{(2)}`$ (the scattered/outgoing field expansion coefficients will differ). In practice, the multipole expansion will be terminated at some $`n=N_{\mathrm{max}}`$. For the case of multipole fields produced by an antenna that is contained within a radius $`a`$, $`N_{\mathrm{max}}=ka`$ is usually adequate, but $`N_{\mathrm{max}}=ka+3\sqrt[3]{ka}`$ is advisable if higher accuracy is needed . This can also be used as a guide for choosing $`N_{\mathrm{max}}`$ for beams—if the beam waist is contained in a radius $`a`$, this can be used to choose $`N_{\mathrm{max}}`$. Since the multipole field necessarily deviates from the paraxial field, the required accuracy is open to question. It appears that $`a=w_0`$ generally gives adequate results.

For focal plane matching of the fields, the regular wavefunctions are used, while for far field matching, the incoming wavefunctions are used. The use of only the incoming portion of the field for matching the fields in the far field means that no knowledge of the outgoing portion of the fields is required.

If the field expansion is truncated at $`n=N_{\mathrm{max}}`$, equation (9) or (10) contain $`2(N_{\mathrm{max}}^2+2N_{\mathrm{max}})`$ unknown variables—the expansion coefficients $`a_{nm}`$ and $`b_{nm}`$. They can be found from a set of $`2(N_{\mathrm{max}}^2+2N_{\mathrm{max}})`$ equations. If we know the field at enough points in space, either by using the standard paraxial form of the beam, or from measurements, a sufficiently large set of equations can be generated. Since the electric field has three vector components at each point in space, three equations result from each point. For robustness, extra points can be used to give an over-determined linear system which can then be solved in a least-squares sense. This is the basic point-matching procedure. The point spacing along any spherical coordinate must be sufficient to prevent aliasing of the VSWFs with respect to that coordinate.

The linear system can be written as a matrix of size $`O(N_{\mathrm{max}}^2)`$, so solution of the system will scale as $`O(N_{\mathrm{max}}^6)`$. While this is clearly less desirable than integral methods which scale as $`O(N_{\mathrm{max}}^4)`$ for surface methods and $`O(N_{\mathrm{max}}^5)`$ for volume methods, integral methods require a much more closely spaced grid of points. Since we are mainly interested in strongly focussed beams, when $`N_{\mathrm{max}}`$ is small, the gain in speed through the use of a coarse grid results in improved performance. We can note that the most efficient solution to very wide beams is simply to approximate them as plane waves (an $`O(N_{\mathrm{max}})`$ solution), which will be quite adequate if the scattering particle is small compared to the beam width, and located near the beam axis. Optimisations for axisymmetric beams, giving $`O(N_{\mathrm{max}}^3)`$ performance, are described later.

## 3 Focal plane matching

Point-matching in the focal plane offers a number of advantages. Firstly, the paraxial beam has a simple phase structure in the focal plane. Secondly, the focal plane is likely to be the region in which scattering or trapping takes place. Thirdly, the irradiance distribution in the focal plane may be known experimentally. Although an accurate measurement of the focal plane irradiance distribution is usually quite difficult, due to the finite resolution of typical detectors, the structure of the beam, and the approximate focal spot size, can be determined. The major disadvantage is that the greatest differences between the paraxial beam and the real beam can be expected to occur in the focal plane. For example, a paraxial beam has zero longitudinal components of the electric field, but, from Maxwell’s equations, the longitudinal component must be non-zero.

For mathematical convenience, we choose a coordinate system so that the beam axis coincides with the $`z`$ axis, with the beam propagating in the $`+z`$ direction, with the focal plane coincident with the $`xy`$ plane (the $`\theta =\pi /2`$ plane). Next, we can generate a grid of points, and calculate the electric field components of the paraxial beam at these points.

A practical difficulty immediately presents itself, however, since the spherical harmonics $`Y_n^m(\theta ,\varphi )`$ with odd $`n+m`$ are odd with respect to $`\theta `$, and will be zero for $`\theta =\pi /2`$—the focal plane. This is to be expected physically, since the same electric field amplitude in the focal plane can correspond to a beam propagating in the $`+z`$ direction, the $`z`$ direction, or a standing wave created by the combination of these. A mathematically rigorous way of dealing with this is to match the derivatives of the electric field amplitude, as well as the field amplitude itself. However, this isn’t necessary, since we have already chosen a direction of propagation as one of the initial assumptions. Therefore, it is sufficient, and computationally efficient, to use only those VSWFs that will be non-zero to determine the corresponding “missing” expansion coefficients. This will give half of the required expansion coefficients. Then, from the direction of propagation, $`a_{nm}=b_{nm}`$ or $`b_{nm}=a_{nm}`$ gives the “missing” coefficients. If the beam propagates in the $`z`$ direction, $`a_{nm}=b_{nm}`$ or $`b_{nm}=a_{nm}`$ are used. If the beam was a pure standing wave combination of counter-propagating beams, these expansion coefficients would be zero.

In principle, this procedure can be carried out for an arbitrary beam. However, determining the vector components of beams where the Poynting vector is not parallel to the $`z`$ axis in the focal plane will be problematic, since there will be no unique choice of electric field direction since the paraxial beam is also a scalar beam. Therefore, we will restrict ourselves to beams with $`z`$-directed focal plane Poynting vectors only. This eliminates, for example, Laguerre-Gaussian beams LGpl with $`l0`$.

Often, the beam of interest will have an axisymmetric irradiance distribution. Any beam that can be written as a sum of $`l=0`$ Laguerre-Gaussian modes LGp0 will satisfy this criterion. The most significant such beam is the TEM00 Gaussian beam (which is also the LG00 beam). Beams with these properties have only $`m=\pm 1`$ multipole components , which follows simply from the $`\mathrm{exp}(\mathrm{i}m\varphi )`$ dependence of the VSWFs. This allows a significant computational optimisation, since only the $`m=\pm 1`$ VSWFs need to be included in the point-matching procedure. This results in $`O(N_{\mathrm{max}}^3)`$ performance.

The point-matching procedure described above was implemented using MATLAB . Since the main computational step is the solution of an overdetermined linear system, a linear algebra oriented language/system such as MATLAB was natural choice. Additionally, existing MATLAB routines for the calculation of Bessel, Hankel, and associated Legendre functions allow simple calculation of VSWFs. Since the VSWFs only need to be calculated for $`\theta =\pi /2`$, a special purpose fast VSWF calculation routine could be used.

### 3.1 TEM00 beam

The scalar amplitude $`U`$ of a TEM00 Gaussian beam is given in cylindrical coordinates by

$$U=U_0\frac{q_0}{w_0+q_z}\mathrm{exp}\left[\mathrm{i}k\left(z+\frac{r^2}{2q_z}\right)\right]$$
(11)
where $`k`$ is the wavenumber, $`w_0`$ is the waist radius, $`U_0`$ is the central amplitude, $`q_z=q_0+z`$ is the complex radius of curvature, and $`q_0=\mathrm{i}z_R`$ where $`z_R=kw_0^2/2`$ is the Rayleigh range. In the focal plane, $`z=0`$ and the above expression can be slightly simplified. For a beam linearly polarised in the $`x`$ direction, $`E_x=U`$ and $`E_y=0`$, while for circularly polarised beams, $`E_x=U`$ and $`E_y=\pm \mathrm{i}U`$. $`E_z=0`$ in all cases.

The general structure of a real beam point-matched with a TEM00 Gaussian beam of waist radius $`w_0=0.5\lambda `$ is shown in figure 1. The beam is very similar to those given by Gouesbet’s localised approximation , although the (very small) $`y`$ component is somewhat different. The robustness of the algorithm is shown by the generation of the correct $`z`$ components for the final beam, even though the original paraxial beam was unphysically assumed to have $`E_z=0`$.

The behaviour of the beam as it is more tightly focussed is shown in figure 2. The increasing axial asymmetry of plane-polarised beams as they are more strongly focussed can be clearly seen. The other feature to note is that focussed beams have a minimum waist size—the more strongly focussed beams shown here closely approach this minimum size, the well-known diffraction limited focal spot . The lack of any such limit in the paraxial beam is blatantly unphysical, and we cannot expect a beam matched to a $`w_0=0.1\lambda `$ paraxial beam to have such an unreasonably small waist radius.

The $`N_{\mathrm{max}}`$ for these calculations was determined by assuming that an effective radius of $`3w_0`$ contains the entire beam in the focal plane. This gives $`N_{\mathrm{max}}=16`$, 9, and 6 for the three beam widths of $`w_0=0.5\lambda `$, $`w_0=0.2\lambda `$, and $`w_0=0.1\lambda `$. All VSWFs were included in the solution (that is, the axisymmetric beam optimisation was not used), so the number of unknown variables for each case was 576, 198, and 96. The fields were matched at a number of points equal to the number of variables, distributed on a polar grid. This is an overdetermined system, with 50% more points than required. The times required for the solution of the linear system on a 1.5 GHz PC were 11.2 s, 0.50 s, and 0.046 s, respectively. If the solution is restricted to the $`m=\pm 1`$ VSWFs, the number of variables is reduced to $`4N_{\mathrm{max}}`$, and the times required for the three cases shown are reduced to 0.006 s, 0.002 s, and 0.001 s.

### 3.2 Bi-Gaussian beam

We consider a simple non-axisymmetric beam, constructed by replacing $`r^2`$ in the paraxial Gaussian beam formula (11) by $`(ax)^2+(by)^2`$, giving a bi-Gaussian irradiance distribution with an elliptical focal spot. The expansion coefficients are no longer restricted to $`m=\pm 1`$.

The irradiance distributions of point-matched bi-Gaussian beams are shown in figure 3. The effect of the minimum spot size is evident, as it first restricts the focussing in the $`y`$ direction, and then the $`x`$ direction, as the beam is more strongly focussed. The odd $`m`$ multipole components are non-zero in these beams. The increased width of the beam results in longer calculations using the $`3w_0`$ cutoff radius. A smaller cutoff radius equal to $`w_0`$ produces very similar results, and may be more suitable for practical use. Since the non-zero expansion coefficients can be predicted from the symmetry of the beam, it would also be possible to speed up the calculations by only including those coefficients.

## 4 Far field matching

In the far field, the beam becomes an inhomogeneous spherical wave. The radial component of the electric and magnetic fields vanishes. The large radius limits for the VSWFs :

$`𝐌_{nm}^{(1,2)}(k\mathrm{r})|_{krn^2}`$ $`=`$ $`{\displaystyle \frac{N_n}{kr}}(\mathrm{i})^{n+1}\mathrm{exp}(\pm \mathrm{i}kr)𝐂_{nm}(\theta ,\varphi )`$ (12)
$`𝐍_{nm}^{(1,2)}(k\mathrm{r})|_{krn^2}`$ $`=`$ $`{\displaystyle \frac{N_n}{kr}}(\mathrm{i})^n\mathrm{exp}(\pm \mathrm{i}kr)𝐁_{nm}(\theta ,\varphi )`$ (13)
can be usefully employed. The far field is best written in spherical polar coordinates, since this gives at most two non-zero components. In cartesian coordinates, all three components will generally be non-zero, although only two will be independent. Since only two vector components of the field can be used in the point matching procedure, the minimum number of points at which to match the fields that is needed is higher.

The non-zero spherical polar vector components of the far field beam can be obtained from the paraxial scalar expressions for the two plane polarised components $`E_x`$ and $`E_y`$:

$`E_\theta `$ $`=`$ $`E_x\mathrm{cos}\varphi E_y\mathrm{sin}\varphi ,`$ (14)
$`E_\varphi `$ $`=`$ $`E_x\mathrm{sin}\varphi +E_y\mathrm{cos}\varphi .`$ (15)
### 4.1 TEM00 beam

The far field spherical wave limit of the paraxial Gaussian beam formula (11) can be found by converting the cylindrical coordinates to spherical coordinates and taking the large $`r`$ limit, giving

$$U=U_0\mathrm{exp}(k^2w_0^2\mathrm{tan}^2\theta /4)$$
(16)
Focal plane irradiance contours for far field matched beams are shown in figure 4. Except for the use of the far field for matching, the procedure, and the results, are the same as the focal plane matching case.

### 4.2 LGpl beams

We recall that the focal plane matching procedure is greatly complicated when the Poynting vector in the focal plane is not parallel to the beam axis, and note that the same does not apply to far field matching—the Poynting vector is always purely radial in the far field regardless of its behaviour in the focal plane. Therefore, beams of these types, such as Laguerre-Gaussian LGpl doughnut beams are best dealt with by far field matching. (The failure of naive attempts to model LGpl beams by focal plane matching simply by using the correct irradiance distribution with the correct $`\mathrm{exp}(\mathrm{i}m\varphi )`$ azimuthal phase variation for $`E_x`$ is readily observed when attempted.) The far field limit for Laguerre-Gaussian beams  is

$$U=U_0(2\psi )^{(l/2)}L_p^l(2\psi )\mathrm{exp}(\psi +il\varphi )$$
(17)
where $`\psi =k^2w_0^2\mathrm{tan}^2\theta /4`$ and $`L_p`$ is the generalised Laguerre polynomial.

### 4.3 Truncation by an aperture

Truncation by an aperture is readily included by restricting non-zero values of the incoming beam to less than the angle of the aperture. If a hard-edged aperture is assumed, a much higher $`N_{\mathrm{max}}`$ is required to accurately truncate the field. The aperture must sufficiently far away for the far field limit to be accurate, which requires a distance $`rn^2/k`$. A soft-edged aperture may also more realistically model the actual physical system. The effect of increasing truncation of a beam is shown in figure 6. The beam is a circularly polarised TEM00 beam of waist radius $`w_0=0.2\lambda `$. Axisymmetry was assumed, with $`N_{\mathrm{max}}=48`$, with the solution of the linear system requiring 0.79 s on a 1.5 GHz PC. The radiation patterns show that the hard-edged aperture is modelled with reasonable, but not perfect accuracy. The errors due to using a finite number of VSWFs to model the sharp edge are exactly analogous to those seen when using a finite number of Fourier terms to model a sharp step. The increase in focal spot size due to diffraction is clearly shown.

## 5 Conclusion

The point-matching method can be successfully used to obtain multipole expansion equivalents of focussed scalar paraxial beams. Multipole expansions are required for scattering calculations using the T-matrix method or GLMT, or calculations of optical forces and torques using these method, and to be able to obtain a satisfactory expansion of a strongly focussed laser beam that is equivalent in a meaningful way to a standard paraxial laser beam is highly desirable. The method is fastest when the beams are strongly focussed, since the maximum degree $`N_{\mathrm{max}}`$ required for convergence is smaller. The method appears usable even when the beam is focussed to the maximum possible extent. Truncation of the beam by apertures is readily taken into account. If it is necessary to truncate the multipole expansion at $`N_{\mathrm{max}}`$ less than the value ideally required for exact representation of the beam (for example, if the T-matrix of the scatterer is truncated at this $`N_{\mathrm{max}}`$), the multipole expansion given by the point-matching method will be the best fit, in a least squares sense, obtainable for this $`N_{\mathrm{max}}`$. If $`N_{\mathrm{max}}`$ is sufficiently large, the multipole expansion will be well-convergent, and the multipole field can be made to be arbitrarily close to the far field/focal plane field with which it is matched. While the multipole beam will not be equal to the original paraxial beam over all space (only on the matching surface), the multipole beam can be considered as a non-paraxial version of the standard beam.

While the emphasis in this paper has been on obtaining multipole expansions of standard beams, the method used is applicable to arbitrary beams, including analytical forms of beams with corrections for non-paraxiality. Since multipole expansions of such beams are still required for the types of scattering calculations considered here, the point-matching method may prove useful applied to such beams. The only restriction on the beam in question is that it must be possible to calculate the fields at a suitable representative set of points.

Compared with integral methods, point-matching is faster for strongly focussed beams since the method tolerates a wider grid point spacing. The chief disadvantage, compared with integral methods, is worse performance for sufficiently large $`N_{\mathrm{max}}`$.

Compared with the localised approximation , point-matching is slower, but is applicable to extremely focussed beams, and allows simple calculation for arbitrary beams, including beams with no known analytical representation.

Practical applications typically require multipole expansions of the beam in a coordinate system centred on a scattering particle, rather than the beam waist as done in the examples presented here. There are two distinct methods in which particle-centred expansions can be determined. Firstly, it should be noted that there is no requirement in the point-matching method for the beam waist to coincide with the $`xy`$ plane and the beam axis to coincide with the $`z`$ axis. Therefore, coordinate axes can be chosen to coincide with the scattering particle, and the points in the focal plane or far field at which the fields are matched specified in this particle-centred coordinate system. This is simply done for focal plane matching. For far field matching, we note that translation of the coordinate system by $`𝐱`$ is equivalent to a phase shift of $`k𝐱\widehat{𝐫}`$ at the points in the far field. A larger $`N_{\mathrm{max}}`$ will typically be required for convergence since a larger radius is required to contain the beam waist when the beam waist is not centred on the origin.

Alternately, and better for repeated calculations, the rotation transformation for VSWFs

$$𝐌_{nm}^{(1,2)}(k𝐫)=\underset{m^{}=n}{\overset{n}{}}D_{m^{}m}^n(\alpha \beta \gamma )𝐌_{nm^{}}^{(1,2)}(k𝐫^{})$$
(18)
and similarly for $`𝐍_{nm}^{(1,2)}(k𝐫)`$ can be used, where $`D_{m^{}m}^n(\alpha \beta \gamma )`$ are Wigner $`D`$ functions , along with the translation addition theorem .

We also note that the paraxial beam waist is a rather misleading parameter to use the describe the multipole beam given by the point-matching procedure, since the actual beam wasit of the multipole beam will differ from that of the paraxial beam. To be able to caclulate a multipole beam of a given waist radius from a paraxial beam, it is necessary to know what paraxial beam waist corresponds to the actual beam waist. Graphs comparing the paraxial and observed waists are given in figure 7. For computational convenience, approximation formulae of the form

$$w_{0\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{a}\mathrm{l}}=w_0+\frac{c_1}{w_0}+\frac{c_2}{w_0^2}+\frac{c_3}{w_0^3}+\mathrm{}$$
(19)
can be used to determine the required paraxial beam waist. The approximation coefficients are listed in table 1.

In conclusion, a reasonable solution to the problem of multipole representation of strongly focussed laser beams is presented. Although the VSWF expansions are not identical to the paraxial beams from which they are derived, they are related in a natural manner. These beams are of particular interest for optical trapping and scattering by single particles within optical traps.