Neutral and non-neutral collisionless plasma equilibria for twisted flux tubes: The Gold-Hoyle model in a background field

We calculate exact one-dimensional collisionless plasma equilibria for a continuum of flux tube models, for which the total magnetic field is made up of the `force-free' Gold-Hoyle magnetic flux tube embedded in a uniform and anti-parallel background magnetic field. For a sufficiently weak background magnetic field, the axial component of the total magnetic field reverses at some finite radius. The presence of the background magnetic field means that the total system is not exactly force-free, but by reducing its magnitude the departure from force-free can be made as small as desired. The distribution function for each species is a function of the three constants of motion; namely the Hamiltonian and the canonical momenta in the axial and azimuthal directions. Poisson's Equation and Amp\`{e}re's Law are solved exactly, and the solution allows either electrically neutral or non-neutral configurations, depending on the values of the bulk ion and electron flows. These equilibria have possible applications in various solar, space and astrophysical contexts, as well as in the laboratory.


I. INTRODUCTION
There has been significant recent work on Vlasov-Maxwell (VM) equilibria that are consistent with nonlinear force-free [1][2][3][4][5][6][7][8] and 'nearly force-free' 9 magnetic fields in Cartesian geometry. Therein, force-free refers to a magnetic field for which the associated current density is exactly parallel, which is the definition we shall also use, These works consider one-dimensional (1D) collisionless current sheets, with Refs. 1-8 specifically calculating VM equilibrium distribution functions (DFs) that are self-consistent with a given specific magnetic field configuration. A natural question to consider is whether it is also possible to find self-consistent force-free (or nearly force-free) VM equilibria for other geometries, in particular cylindrical geometry. In this paper we shall present particular VM equilibria for 1D magnetic fields which are nearly force-free in cylindrical geometry, i.e. flux tubes/ropes Two of the archetypal field configurations in cylindrical geometry are the z-Pinch (with axial current and azimuthal magnetic field), a classical example of which is the Bennett Pinch 10 ; and the θ-Pinch (azimuthal current and axial magnetic field). Consideration of 'Vlasovfluid' models of z-Pinch equilibria was given in Ref. 11, with Ref. 12 calculating z-Pinch equilibria and an extension with azimuthal ion-currents. Others have also constructed kinetic models of the θ-pinch, see Refs. 13 and 14 for examples. In the same year as Pfirsch 15 , cylindrical kinetic equilibria with only azimuthal current were studa) Electronic mail: oliver.allanson@st-andrews.ac.uk ied in Ref. 16 . For examples of treatments of the stability of fluid and kinetic linear pinches, see Refs 15,17, and 18 respectively.
Recently there have been studies on 'tokamak-like' VM equilibria with flows [19][20][21] , starting from the VM equation in cylindrical geometry and working towards Grad-Shafranov equations for the vector potential. We also note two Vlasov equilibrium DFs in the literature that are close in style to the one that we shall present. The first is described in a brief paper 22 , with an equilibrium presented for a cylindrical pinch. However, their distribution describes a different magnetic field and the DF appears not to be positive over all phase space. The second DF is a very recent paper 23 that actually describes a magnetic field much like the one that we discuss. Their DF is designed to model 'ion-scale' flux tubes in the Earth's magnetosphere. Formally, their quasineutral model approaches a nonlinear force-free configuration in the limit of a vanishing electron to ion mass ratio. In their model, current is carried exclusively by electrons and the nonnegativity of the DF depends on a suitable choice of microscopic parameters. Finally, we mention that in beam physics, much work on constructing cylindrical VM equilibria is done by looking for mono-energetic distributions with conserved angular momentum, see Refs. [24][25][26][27] for some examples.
Magnetic flux tubes and flux ropes are prevelant in the study of plasmas, with a wide variety of observed forms in nature and experiment, as well as uses and applications in numerical experiments and theory. Some examples of the environments and fields of study in which they feature include solar 28,29 ; solar wind 30,31 ; planetary magnetospheres 32,33 and magnetopauses 34 ; astrophysical plasmas 35,36 ; tokamak 37,38 , laboratory pinch experiments 39 and the basic study of energy release in magnetised plasmas 40 , to give a small selection of references.
One application of flux tubes is in the study of solar active regions 41 and the onset of solar flares and coronal mass ejections [42][43][44] . A classic magnetohydrodynamic (MHD) model for magnetic flux tubes was first presented by T. Gold and F. Hoyle (GH) 45 , initially intended for use in the study of solar flares. The GH model is an infinite, straight, 1D and nonlinear force-free magnetic flux tube with constant 'twist' 46 . Mathematically, the GH magnetic field could be regarded as the cylindrical analogue 47 of the Force-Free Harris sheet 2 (a planar current sheet model), as the Bennett Pinch 10 might be to the 'original' Harris Sheet 48 . It is typical to consider solar, space and astrophysical flux tubes within the framework of MHD, e.g. see Ref. 49. However, many of these plasmas can be weakly collisional or collisionless, with values of the collisional free path large against any fluid scale 50 , making a description using collisionless kinetic theory necessary. It is the intention of this paper to study the GH flux tube model beyond the MHD description, since -apart from the very recent work in Ref. 23 -we see no attempt in the literature of a microscopic description of the GH field. Other than any interesting theoretical advances, a possible application of the results of this study could be to implement the obtained model in kinetic (particle) numerical simulations.
In Cartesian geometry, the work in Refs. 1-8 used the method proposed by Channell 51 to tackle the VM inverse problem, i.e. to determine self-consistent equilibrium DFs for a given magnetic field configuration. Channell described the extension of his work to cylindrical geometry as 'not possible in a straightforward manner.' As explained in Ref. 20 (in which cylindrical coordinates are used to model a torus), this is due in part to the 'toroidicity' of the problem, i.e. the 1/r factor in the equations. As we shall see in this paper, another potential complication is the need to allow -at least in principle -a non-zero charge density. The work in this paper does not present a generalised method for the VM inverse problem in cylindrical geometry, but instead some particular solutions for a specific given magnetic field.
The paper is structured as follows. In Section 2 we first review the theory of the equation of motion consistent with a collisionless DF in cylindrical geometry, and discuss the question of the possibility of 1D force-free equilibria. Then we introduce the magnetic field to be used. We note that whilst the work in this paper is applied to a particular magnetic field from Subsection 2.1 onwards, the steps taken to calculate the equilibrium DF seem as though they could be adaptable to other cases. In Section 3 we present the form of the DF that gives the required macroscopic equilibrium, and proceed to 'fix' the parameters of the DF by explicitly solving Ampère's Law and Poisson's Equation. Note that whilst we choose to consider a two-species plasma of ions and electrons, we see no obvious reason preventing the work in this paper being used to describe plasmas with a different composition. In Section 4. we present a preliminary analysis of the physical properties of the equilibrium. Particularly technical calculations are in the Appendices. Appendix A contains the zeroth and first order moment calculations, used to find the number densities and bulk flows directly, and in turn the charge and current densities. Appendix B contains the mathematical details of the existence and location of multiple maxima of the DF in velocity-space.

A. The Vlasov equation and the equation of motion
A collisionless equilibrium is characterised by the 1particle distribution function, f s , a solution of the steadystate Vlasov equation (e.g. see Ref. 52). The Vlasov equation in cylindrical coordinates is see for example Refs. 16, 19, and 53. Here i, j and k are used as 'spatial' indices running over {1, 2, 3}, and s is used as the particle species index. Individual particle positions and velocities are given by ( respectively, for r the horizontal distance from the z axis, and θ the azimuthal angle. The totally antisymmetric unit tensor of rank 3 (the Levi-Civita tensor) is ε ijk , and the Einstein summation convention is applied (such that repeated indices are summed over, with subscript and superscript indices used to describe co-and contravariant components respectively. The mass and charge of particle species s are m s and q s respectively. The electric and magnetic fields are defined as E = −∇φ and B = ∇ × A, for φ the scalar potential and vector potential A. The 'fluid' equation of motion of a particular species s is found by taking first-order velocity moments of the Vlasov equation. After a routine but laborious momenttaking calculation, we see that -in equilibrium (∂/∂t = 0), assuming a one-dimensional configuration with only radial dependence (∂/∂θ = ∂/∂z = 0), and letting f s be an even function of the radial velocity v r -force balance for species s is maintained according to The pressure tensor for species s is a rank-2 tensor and is defined by where v i = u is + w is , for u is the bulk velocity of species s and v i the individual particle velocity. Note that the assumption of f s to be an even function of v r automatically implies that u rs = P rθ = P zr = 0. Equation (2) can be summed over species to give where is the force density associated with the rotating bulk flows of the ions and electrons. Equation (3) is a cylindrical analogue of the force balance equation in Cartesian geometry (e.g. see 54). There are 'extra inertial terms' as compared to the case of Cartesian geometry. From the point of view of a particular magnetic field B (which is the point we take by specifying a particular macroscopic equilibrium), we see that equilibrium is maintained by a combination of density/pressure variations as in the case of Cartesian geometry, but with additional contributions from centrifugal forces and as an inevitable result of the resultant charge separation, an electric field. This clearly demonstrates that 'sourcing' an exactly force-free macroscopic equilibrium with an equilibrium DF in a 1D cylindrical geometry is inherently a more difficult task than in the Cartesian case. The presence of 'extra' positive definite inertial forces and, almost inevitably, forces associated with charge separation raises the question of whether exactly force-free equilibria are possible at all in this paradigm. Before proceeding, we comment that given certain macroscopic constraints on the electromagnetic fields or fluid quantities -such as the force-free condition, or a specific given magnetic field (for example) -it is not a priori known how to calculate a self-consistent Vlasov equilibrium, or if one even exists within the framework of the assumptions made. Hence one has to proceed more or less on a case by case basis, with the intention of achieving consistency with the required macroscopic conditions, upon taking moments of the DF.

B. Methods for calculating an equilibrium DF
In Refs. 2 and 51 for example, a method used to calculate a DF, given a prescribed 1D magnetic field was Inverse Fourier Transforms (IFT). A distribution function of the form f s ∝ e −βsHs g s (p xs , p ys ), was used, with H s , p xs and p ys the conserved particle Hamiltonian and canonical momenta in the x and y directions, and g s an unknown function, to be determined. Since our problem is one of a 1D equilibrium with variation in the radial direction, the three constants of motion are the Hamiltonian, and the canonical momenta in the θ and z directions: A function of a subset of the constants of motion is automatically a solution of the VM equation (e.g. see Ref. 52). One can try to calculate an equilibrium distribution for the Gold-Hoyle force-free flux tube without a background field by a similar method, assuming a DF of the form f s ∝ e −βsHs g s (p θs , p zs ).
By exploiting the convolution in the definition of the current density, Ampère's law can be solved by IFT, with the quantity p s defined by p rs = p rs , p θs = p θs r , p zs = p zs .
Notice how when written in this integral form, j is not only a function of A, but -in contrast with the Cartesian case -also of the relevant spatial co-ordinate, r. In the case of zero scalar potential, the result of the calculation is to give a distribution function that is not a solution of the Vlasov equation as it is not a function of the constants of motion only. In essence, an additional exp(−r 2 ) factor is required in the DF to counter exp(r 2 ) terms that manifest by completing the square in the integration. The physical cause here would appear to be the inertial forces associated with the rotational bulk flow. If one assumes a non-zero scalar potential, then it seems impossible to satisfy Ampère's Law. The physical cause seems to be that, in the case of force-free fields, one would require a 'different' electrostatic potential to balance the inertial forces for the ions and electrons, which is of course nonsensical. Thus, our investigation seems to suggest that it is not possible to calculate a DF of the form of equation (6) for the exact GH field.

C. The magnetic field: A Gold-Hoyle flux tube plus a background field
To make progress, we introduce a background field in the negative z direction. The mathematical motivation for this change is to balance the 'exp(r 2 ) problem'. Physically, it seems that the background field introduces an extra term (whose sign depends on species) into the forcebalance, to allow for both the ion and electrons to be in force balance simultaneously, given one unique expression for the scalar potential.
The vector potential, magnetic field and current den-sity used in this paper are as follows (GH+B): B(r) = B 0 0,r 1 +r 2 , The magnetic permeability in vacuo is given by µ 0 and the characteristic magnetic field strength by B 0 . The constant τ has units of inverse length, and we use 1/τ to represent the characteristic length scale of the system (r = τ r) (see Table 1 for a concise list of the dimensionless quantities used in this paper, all denoted with a tilde,˜). The dimensionless constant k > 0 controls the strength of the background field in the z direction, and as a result there are now two different interpretations to be made. We could either consider the system as a GH flux tube of uniform twist embedded in an untwisted uniform background field, or consider the whole GH+B magnetic field as a non-uniformly twisted flux tube. We note that flux tubes embedded in an axially directed background field have recently been observed during reconnection events in the Earth's magnetotail, by the Cluster spacecraft 55 .
In the first interpretation, τ is a direct measure of the 'twist' of the embedded flux tube (see Ref. 46), with the number of turns per unit length (in z) along a field line given by τ /(2π) 45 . In the second interpretation, we see that the system is not uniformly twisted, with the z distance traversed when following a field line (e.g. 56) given by The fact that this depends on r demonstrates that the system as a whole has non-uniform twist. The number of turns per unit length in z of the GH+B field: the 'twist' is given by and is plotted in Figure 1 for three values of k. Since k < 1/2 corresponds to the field-reversal regime, we see a mixture of positive and negative twists ( Figure 1a). However, for k ≥ 1/2 we see only negative values of the twist (Figures 1b and 1c), i.e. we travel in the negative z direction as we wind round the GH+B flux tube in the anti-clockwise direction. The magnetic field is plotted in Figures 2a-2b for two values of k. The k = 0.3 case contains a reversal of theB z field direction and as such is akin to a Reversed Field Pinch (e.g. see Ref. 57 for a laboratory interpretation): this configuration may be of use in the study of astrophysical jets, see Ref. 36 for example. The value k = 1/2 gives zeroB z atr = 0, and as such is the value that distinguishes the two different classes of field configuration, namely unidirectional (k ≥ 1/2) or including field reversal (k < 1/2). The value ofr for which theB z field reverses is plotted in Figure 2c. The magnitude of the GH+B magnetic field is plotted in Figure 3 for three values of k. For all values of k, |B| → 2k for larger, i.e. to a potential field.
The primary task of this paper is to calculate selfconsistent collisionless equilibrium distribution functions for the GH+B field. This problem essentially reduces to solving Ampère's Law such that equation (1) is satisfied. We assume nothing about the electric field however, and in fact use that degree of freedom to solve Ampère's Law. The resultant form of the scalar potential is then substituted into Poisson's equation, to establish the final relationships between the microscopic and macroscopic parameters of the equilibrium.

III. THE EQUILIBRIUM DISTRIBUTION FUNCTION
Although the IFT method did not yield a selfconsistent equilibrium DF for the GH field without a background field, the outcome of the calculation can still be used as an indication of possible forms for the DF for the GH+B field. Using trial and error we arrived at the distribution function which is a superposition of two terms that are consistent macroscopically with a 'Rigid-Rotor', see Ref. 18 for example. A Rigid-Rotor is microscopically described by a DF of the form F (H − ωp θ − V p z ) (with V = 0 in the second term of the DF in equation (10)). Each F (H −ωp θ −V p z ) term corresponds to an average macroscopic motion of rigid rotation with angular frequency ω, and rectilinear motion with velocity V . The dimensionless constantsω s ,Ũ zs ,Ṽ zs and C s are yet to be determined, with C s > 0 for positivity of the distribution. Note that the thermal beta is β s = 1/(k B T s ) and v th,s is the thermal velocity of species s. The ratio of the thermal Larmor radius, r L = m s v th,s /(e|B|) (for e = |q s |) to the macroscopic length scale of the system L(= 1/τ ), is given by typically known as the 'magnetisation parameter' 58 . In our system, the magnitude of the magnetic field and hence δ s itself is spatially variable. For the purposes of the calculations in this paper however, we set as a characteristic value (see Table II for a concise list of the micro and macroscopic parameters of the equilibrium).
A. Maxwell's equations: Fixing the parameters of the DF By insisting on a specific magnetic field configuration (the GH+B field) we have made a statement on the macroscopic physics. In searching for the equilibrium DF, we are trying to understand the microscopic physics. In this sense we are tackling an 'inverse problem'. Once an assumption on the form of the DF is made thenshould the assumed form be able to reproduce the correct moments -this inverse problem reduces to establishing the relationships between the microscopic and macroscopic parameters of the equilibrium. In this section we 'fix' the free parameters of the DF in equation (10), such that Maxwell's equations are satisfied; Note that the solenoidal constraint and Faraday's law are automatically satisfied for the GH+B field in equilibrium,

Ampère's Law
In Appendix A we have calculated the j z current density, found by summing first order moments in v z of the DF. We now substitute in the macroscopic expressions for j z (r), A θ (r) and A z (r) from (9) and (7) into the expression for the j z current density of equation (A4). After this substitution, we can calculate a φ(r) that makes the system consistent. The substitution of the known expressions for j z , A z and A θ gives zs +r 2ω2 s )/2−sgn(qs)ωsr 2 k/δs 1 +r 2 sgn(qs)(ωs−Ũzs)/(2δs) +Ṽ zs C s eṼ 2 zs /2 1 +r 2 −sgn(qs)Ṽzs/(2δs) (13) In order to satisfy the above equality we can construct a solution by introducing a 'separation constant' γ 1 = 0, 1. We multiply the above equation by (1 +r 2 ) 2 which makes the left-hand side constant, whilst the right-hand side is a sum of two terms, one depending on ion parameters and the second depending on electron parameters. Then we can define γ 1 by associating the 'ion term' with the first term on the right-hand side of (14), and the 'electron term' with the second term on the right-hand side of (14). After some algebra we can rearrange these two associations to give two expressions for the scalar potential, one in terms of the ion parameters, and one in terms of the electron parameters: ze +r 2ω2 e )/2+ωer 2 k/δe 1 +r 2 2−(ωe−Ũze)/(2δe) +Ṽ ze C e eṼ 2 ze /2 1 +r 2 2+Ṽze/(2δe) The two values of the scalar potential above must be made identical by a suitable choice of relationships between the ion and electron parameters. Given enough freedom in parameter space, we could say that the z component of Ampère's Law is implicitly solved the above equations, in that one just needs to choose a consistent set of parameters. However, we seek a solution in an explicit sense.
In order to make progress we non-dimensionalise the above equations by multiplying both sides by eβ r with Once this is done we can write the scalar potential in the form Specifically, equations (15) and (16) require the equality of the arguments of the logarithm to hold in order for a meaningful solution to be obtained for the scalar potential. A first step towards this is made by requiring consistent powers of the 1 +r 2 'profile' in the right-hand side of the above expression to allow factorisation. Hence and hence the rigid-rotation,ω s , is fixed by the difference of the rectilinear motion,Ũ zs −Ṽ zs . On top of this, we require that the power of the 1 +r 2 'profile' on the right-hand side is the same for both the ions and electrons, thus This condition seems to be a statement on an average potential energy associated with the particles. Once more to allow factorisation of the 1 +r 2 'profile', we insist that net exp(r 2 ) terms cancel, i.e.
The physical meaning of this condition seems to be that the frequencies of the rigid rotor for each species are matched according to the relevant magnetisation, and the background field magnitude. The remaining task is to ensure equality of the 'coefficients' These seem to be conditions on the ratios of the energy densities associated with the bulk rectilinear motion and the magnetic field respectively. Thus far we have 8 constraints and 12 unknowns (Ũ zs ,Ṽ zs ,ω s , C s , n 0s , β s ) given fixed characteristic macroscopic parameters of the equilibrium B 0 , τ , and k. We can now write down an expression for φ that explicitly solves the z component of Ampère's law; with φ(0) = 1 eβ r ln D.
Clearly, we require that D > 0 for the expression above to make sense. It is clear that the sign of γ 1 could, in principle, affect the sign of D. It is seen from (20) that positivity of D implies that By rearranging the above inequalities to make C s the subject, it can be seen after some algebra that positivity of D and C s is guaranteed when Note that these conditions are sufficient, but not necessary, i.e. it is possible to have D > 0 and C s > 0 for any value of γ 1 = 0, 1, and even for sgn(Ũ zs ) = sgn(Ṽ zs ) in the case of γ 1 < 0. Thus far we have only considered the j z component, and it is premature to consider all components of Ampère's Law satisfied. Let us move on to consider the θ component. In a process similar to that above, we substitute in the macroscopic expressions for j θ (r), A θ (r) and A z (r) for the GH+B field into the expression for the j θ current density of equation (A6) in Appendix A. After this substitution, we can once more calculate the φ that makes the system consistent. The substitution gives Using the parameter relations as above, we determine that the scalar potential is again given in the form of (21), φ(r) = 1 eβ r E ln 1 +r 2 + φ(0).
Hence, this form of the scalar potential is consistent provided for γ 2 = 1 another separation constant. These seem to be conditions on the ratios of the energy densities associated with the bulk rotation and the magnetic field respectively. This has added two more constraints. Once again we must ensure that D > 0. Since ω e < 0, the right-hand side of the above equation implies that γ 2 > 0 to ensure that D > 0. Whilst the left-hand side implies that γ 2 < 1 for positivity of D since ω i > 0. Hence we can say that for positivity 0 < γ 2 < 1.
We can now consider Ampère's Law satsified, given a φ that solves Poisson's equation. As a result, the problem of consistency is now shifted to solving Poisson's Equation, where the remaining degrees of freedom lie.

Poisson's Equation
The electric field is calculated as E = −∇φ, giving We can now take the divergence of the electric field ∇ · E = τr −1 ∂r(rE r ) and so This gives a non-zero net charge per unit length (in z) of The charge density derived in equation (27) The second equality is found by substituting the form of the vector potential from equation (7), and the final equality is reached by using the conditions derived in equations (17) - (21).
We can now match equations (27) and (29) We now have 12 physical parameters (Ũ zs ,Ṽ zs ,ω s , C s , n 0s , β s ) with 11 constraints (17)(18)(19)(20), (25) & (30). For example, if one picks B 0 , τ , k and one microscopic parameter, say β i , then the remaining parameters of the equilibrium, (Ũ zs ,Ṽ zs ,ω s , C s , n 0s , β e ), are now determined. One could of course choose the values of a different set of parameters, and determine those that remain by using the constraints derived. Note that whilst the constants γ 1 = 0, 1 and 0 < γ 2 < 1 are system parameters, they are not physically meaningful as they only represent a change in the gauge of the scalar potential.

A. Non-neutrality & the electric field
It is seen from equations (27) and (28) that basic electrostatic properties of the equilibrium described by f s are encoded in E. The equilibrium is electrically neutral only when E = 0, and non-neutral otherwise. Specifically, there is net negative charge when E > 0, and net positive charge when E < 0. This net charge is finite in the (r, θ) plane and given by Q in equation (28).
Physically, the sign of E seems to be related to the respective magnitudes of the bulk rotation frequencies, ω s . From equations (17) and (18) we see that E > 0 implies thatω and E < 0 implies that Hence, E > 0 is seen to occur for 'sufficiently large' bulk ion rotation frequencies, and 'sufficiently small' (in magnitude) bulk electron rotation frequencies. A positive E corresponds to an electric field directed radially 'inwards'. This seems to make sense physically, by the following argument. A 'larger' (ω i > ω i ) bulk ion rotation freqency gives a 'larger' centrifugal force, and a 'smaller' (|ω e | < ω e ) bulk electron rotation frequency gives a 'smaller' centrifugal force. For a dynamic interpretation, at a fixed r, the ions are forced to a slightly larger radius than the electrons, i.e. a charge separation manifests on small scales. This charge separation results in an inward electric field, E r < 0. An equally valid interpretation is to say that for an equilibrium to exist, an electric field must exist to counteract the differences in the centrifugal forces associated with the bulk ion and electron rotational flows.
In a similar manner, E < 0 is seen to occur for 'sufficiently small' (ω i < ω i ) bulk ion rotation frequencies, and 'sufficiently large' (|ω e | > ω e ) bulk electron rotation frequencies. A negative E corresponds to an electric field directed radially 'outwards'. We can then interpret these result physically, in a manner like that above.
Finally, we can interpret the neutral case, E = 0, as the intermediary between the two circumstances considered above. That is to say that the equilibrium is neutral when the bulk rotation flows are just matched accordingly, such that there is no charge separation and hence no electric field.

B. The equation of state and the plasma beta
For certain considerations, e.g. the solar corona, it would be advantageous if the DF had the capacity to describe plasmas with sub-unity values of the plasma beta: the ratio of the thermal energy density to the magnetic energy density For our configuration, the number density is seen to be proportional to the rr component of the pressure tensor, P rr,s = n s k B T s . This is demonstrated by the following calculation. In order to calculate P rr , we must consider the integral However, we do not have to consider a bulk velocity in the r direction here (u rs = 0), since f s is an even function of v r . Using the fact that and by consideration of equations (32) and the number density, we see that that is to say that k B T s = m s v 2 th,s . Note that if n i = n e := n and hence E = 0 (neutrality), then we have an equation of state given by This resembles expressions found in the Cartesian case, in Refs. 3, 7, and 51 for example. Incidentally, we can use the connection between n s and P rr to give an expression for the β pl that is perhaps more typically seen, The square magnitude of the magnetic field (equation (8)) is given by Using the number density from equation (A1) in the definition of the plasma beta from equation (31), as well as the equilibrium conditions (17) - (21) gives It is not immediately obvious from the above equation what values β pl can have. However it is readily seen that asr → ∞ then β pl → 0, essentially since the number density is vanishing at large radii. On the central axis of the tube we see that suggesting that for a suitable choice of parameters, it should be possible to attain any value of β pl on the axis.

C. Plots of the DF
A characteristic that one immediately looks for in a new DF is the existence of multiple maxima in velocity space, which are a direct indication of nonthermalisation, relevant for the existence of microinstabilities (e.g. see 59 ). Using an analysis very similar to that in 3 , we can derive -for a given value ofω s -conditions onr and eitherṽ z orṽ θ , for the existence of multiple maxima in theṽ θ orṽ z direction respectively. We present these calculations in Appendices B 1 and B 2. The most readily understood results are that multiple maxima in theṽ θ direction can only occur forr > 2/|ω s |, and in thẽ v z direction for |ω s | > 2. Given these necessary conditions, one can then calculate that multiple maxima of f s will occur in theṽ θ direction forṽ z bounded above and below, and vice versa.
In Figures (4-7) we present plots of the DFs over a range of parameter values. Figures (4) and (5) show the ion DFs for k = 0.1 and k = 1 respectively, for all combinations ofω i = 1, 3,r = 0.5, 2 and C s = 0.1, 1, and with the magnetisation parameter δ i = 1. As a graphical confirmation of the above discussion, we can only see multiple maxima in theṽ θ direction forr > 2/|ω s |, and in theṽ z direction for |ω s | > 2, with the appropriate bounds marked by the horizontal/vertical white lines.
Aside from multiple maxima in the orthogonal directions, the DF can also be 'two-peaked'. That is, the DF can have two isolated peaks in (ṽ z ,ṽ θ ) space. This is seen to occur for figures (5d, 5g, 5h). Hence, f i is seen to be 'two-peaked' when k = 1 for bothr > 2/ω i and r < 2/ω i . However, we do not see a two-peaked DF for k = 0.1. This seems to suggest that the stronger guide field (k = 1) correlates with multiple peaks. Physically, this may correspond to the fact that a homogeneous guide field is consistent with a Maxwellian DF centred on the origin in (ṽ z ,ṽ θ ) space, given that a Maxwellian contributes zero current. Hence, if the 'main' part/peak of the DF is centred away from the origin, then the Maxwellian contribution from the guide field could contribute a secondary peak. These secondary peaks are seen to be more pronounced whenC i is larger, i.e. the contribution from the second term from the DF is greater.
Figures (6) and (7) show the electron DFs for k = 0.1 and k = 1 respectively, for all combinations ofω e = 1, 3; r = 0.5, 2, and C e = 0.1, 1, and with the magnetistaion parameter δ e = δ i m e /m i ≈ 1/ √ 1836. This choice of magnetisation corresponds to T i = T e . In general we see DFs with fewer multiple maxima in velocity space than the ion plots, which is physically consistent with the electrons being more magnetised, i.e. more 'fluidlike'. In particular we see no multiple maxima in figure  7, the case with the stronger background field.
Note that when the electrons to have the same magnetisation as the ions, i.e. δ e = δ i = 1, then these marked differences in the velocity-space plots disappear, and we observe a qualitative symmetry f i (ṽ θ ,ṽ z , r) ∝ f e (−ṽ θ , −ṽ z , r).

V. SUMMARY
In this paper we have calculated one-dimensional collisionless equilibria for a continuum of magnetic field models based on the Gold-Hoyle flux tube, with an additional constant background field in the axial direction. This study was motivated by a desire to extend the existing methods for solutions of the 'inverse problem in Vlasov equilibria' in Cartesian geometry, to cylindrical geometry. Initial efforts focussed on solving for the exact forcefree Gold-Hoyle field, but this seems impossible due to the positive definite centrifugal forces. The Gold-Hoyle field in particular was chosen as it represents the 'natural' analogue of the Force-Free Harris Sheet in cylindrical geometry, a magnetic field whose VM equilibria have been the subject of recent study, 2-4,6-8 .
A background field was introduced, and an equilibrium distribution function was found that reproduces the required magnetic field, i.e. solves Ampère's Law. It is the presence of the background field that allows us to solve Vlasov's equation and Ampère's Law, and it appears physically necessary as it introduces an 'asymmetry'; namely an extra term into the equation of motion whose sign depends explicitly on species. In contrast to the 'demands' of insisting on a particular magnetic field, no condition was made on the electric field. The distribution function allows both electrically neutral and nonneutral configurations, and in the case of non-neutrality we find an exact and explicit solution to Poisson's equation for an electric field that decays like 1/r far from the axis. We note here that the type of solutions derived in this paper could -after a Galilean transformation -be interpreted as 1D BGK modes with finite magnetic field (see Refs. 60-63 for example, to provide some context).
An analysis of the physical properties of the DF was given in Section IV, with some detailed calculations in Appendix B. The dependence of the sign of the charge density (and hence the electric field) on the bulk ion and electron rotational flows was analysed, with a physical interpretation given. Essentially the argument states that the electric field exists in order to balance the difference in the centrifugal forces between the two species. The DF was found to be able to give sub-unity values of the plasma beta, should this be required/desirable given the relevant physical system that it is intended to model. The final part of the analysis focussed on plotting the DF in velocity space, for certain parameter values, and at different radii. Mathematical conditions were found that determine whether or not the DF could have multiple maxima in the orthogonal directions in velocity space, and these are corroborated by the plots of the distribution functions. For certain parameter values, the DF was also seen to have two separate, isolated peaks. This non-thermalisation suggests the existence of microinstabilities, for a certain choice of parameters.
Further work could involve a deeper anlysis of the properties of the distribution functions and their stability. This work has also raised a fundamental question: 'is it possible to describe a one-dimensional force-free collisionless equilibrium in cylindrical geometry?' Preliminary investigations seem to suggest that it is not possible. In this appendix we calculate the zeroth and first order velocity space moments of the DF, necessary for the charge density and the current density respectively. See Table 1 for a clarification of all dimensionless quantities denoted by a tilde,˜. The number density of species s is given by the zeroth moment of the DF; Thep rs dependence of the DF is irrelevant to our discussion, and as such can be integrated out. We can also neglect the scalar potential φ. The reduced DF,F s , in dimensionless form is which then reads We have writtenF s in terms of the canonical momenta, and so we search for stationary points given by ∂F s /∂p θs = 0, equivalent to ∂F s /∂ṽ θs = 0. Setting ∂F s /∂p θs = 0 gives p θs −rÃ θs =ω To derive a necessary condition for multiple maxima, we analyse the RHS of equation (B2), R(p θs ). This function is bounded between 0 and A, and is monotonically increasing. Hence, using techniques similar to those in 3 , a necessary condition for multiple maxima in the DF is that This condition can be shown to be equivalent to Aω s /4 > 1 and soω This demonstrates that for sufficiently smallr, there cannot exist multiple maxima. Equivalently, this condition will always be satisfied for somer, and as such is just a condition on the domain, inr, for which multiple maxima can occur. This condition is not sufficient however, as it could still be the case that there exists only one point of intersection (and hence one maximum), depending on the value of B. It is seen that R has unit slope at Clearly R has unit slope for two values ofp θs . After some graphical consideration of the problem, it becomes apparent that B should be bounded above and below for multiple maxima. After elementary consideration of the functional form of (B2), for example with graph plotting software, we see that multiple maxima in theṽ θ direction can only occur, for a givenr, when B (and henceṽ z ) satisfies these inequalities for ions and these for electrons

Maxima of the DF in vz space
We shall once again use the reduced DF defined in equation (B1) in our analysis. Thus, we shall consider ∂F s /∂p zs = 0, which is equivalent to ∂F s /∂ṽ zs = 0. Setting ∂F s /∂p zs = 0 gives p zs −Ã zs =Ũ zs + C sṼzs e −ωs(pzs+p θs ) 1 + C s e −ωs(pzs+p θs ) To derive a necessary condition for multiple maxima, we analyse the RHS of equation (B8). Each R function is bounded and monotonic. Once again using techniques similar to those in 3 , a necessary condition for multiple maxima in the DF is that max pzs (R 1 (p zs ) + R 2 (p zs )) > 1. (B8) After some algebra this condition can be shown to be equivalent toω 2 s /4 > 1 and so This condition is not sufficient however, as it could still be the case that there exists only one point of intersection, depending on the value of B 1 (= 1/B 2 ). The transition between 3 points of intersection and one occurs at the value of B 1 for which the straight line of slope unity throughp zs = 0 just touches R 1 (p zs ) + R 2 (p zs ) at the point where it also has unit slope. It is readily seen that R 1 + R 2 has unit slope at We see again that R has unit slope for two values ofp zs . Once again, after some graphical consideration of the problem, it becomes apparent that B 1 should be bounded above and below for multiple maxima. After elementary consideration of the functional form of (B8), for example with graph plotting software we see that multiple maxima in theṽ z direction can only occur, for a givenr, when B 1 (and henceṽ θ ) satisfies these inequalities for ions and these for electrons  Characteristic magnetic field strength ms Mass of particle τ Measure of the twist of flux tube qs, q Charge, magnitude of charge k Strength of the background field βs = 1/(kBTs) Thermal beta γ1 = 0, 1, 0 < γ2 < 1 Gauge for scalar potential v th,s Thermal velocity Uzs, Vzs Bulk rectilinear flows δs(r), δs

Magnetisation parameters ωs
Bulk angular frequency n0s Normalistaion of particle number   2: 2a and 2b show the GH+B magnetic field in the xy plane, for two values of k. The curved arrows indicate the direction of theB θ components, whilst the blue-black-red shading denotes the magnitude and direction of thẽ B z component. The k = 0.3 case contains a reversal of theB z field direction and as such is a Reversed Field Pinch whilst k = 0.5 gives zeroB z atr = 0. 2c shows the radius at whichB z changes its direction, for 0 < k < 1/2.B z does not reverse for k ≥ 1/2.    Contour plots of f e in (ṽ z ,ṽ θ ) space for an equilibrium without field reversal (k = 1 > 0.5), for a variety of parameters (ω e ,r, C e ) and δ e ≈ 1/ √ 1836. Note that there are not any multiple maxima in this case.