Some general solutions for linear Bragg-Hawthorne equation

Linear cases of Bragg-Hawthorne equation for steady axisymmetric incompressible ideal flows are systematically discussed. The equation is converted to a more convenient form in a spherical coordinate system. A new vorticity decomposition is derived. General solutions for 16 linear cases of the equation are obtained. These solutions can be specified to gain new analytical vortex flows, as examples in the paper demonstrate. A lot of well_known solutions like potential flow past a sphere, Hill's vortex with and without swirl, are included and extended in these solutions. Special relations between some vortex flows are also revealed when exploring or comparing related solutions.


I. INTRODUCTION
Bragg-Hawthorne equation 1 (or "B-H equation" in short) plays a central role in the study of axisymmetric steady flow of incompressible ideal fluids. The solutions of this equation (and its equivalent equation in magnetohydrodynamics, i.e., the Grad-Shafranov equation) can be applied to support a variety of research from hydrodynamics to magnetohydrodynamics and plasma physics, from vortex and turbulence study to rocket engine 2 and tokamak fusion reactor research, 3 from tornado study 4 to astrophysics, 5,6 etc.
Due to complexity of the equation, however, analytically solvable cases are rare. With special assumptions or restrictions for Bernoulli function and the azimuthal velocity, the equation reduces to a linear one. Most of the known analytic solutions so far are for such linear cases.
In this paper, an attempt is made to explore the linear cases in a systematic way and gain new explicit analytical solutions. First, the equation is rewritten in spherical coordinate system and the variable of polar angle is changed from h to cos h. This gives a new form of the equation that is more convenient to solve in most cases. Following that, a series of assumptions for Bernoulli function and the azimuthal velocity that make the equation linear are listed. These assumptions lead to 16 combinations, each associated with a special linear case of the equation. These 16 cases are then solved, mostly by separate variable method. The results are a series of general solutions, which, when specified, bring new analytical solutions of vortex flows. Many wellknown solutions, including the potential flow around a sphere, Hill's vortex with and without swirl, 7, 8 Bogoyavlenskij's solution of Beltrami flow, 9 the Fraenkel-Norbury family of vortex rings, 10,11 etc., can be obtained from these general solutions as well when constants are set to certain values.
Some linear cases of the equation have a close relation, so do their solutions and the related flows. For example, Hill's vortex without swirl is actually a combination of a uniform flow and an extra velocity. Hill's vortex with swirl can be considered as the Beltrami spherical vortex (the eigenvector field of curl operator in the spherical coordinate system, as discussed in Sec. VII) plus a special rotation. These are discussed in detail with their stream functions and velocities.
A new decomposition of vorticity is also derived. This decomposition applies to all flows (that follow B-H equation) and helps us to understand the impacts of Bernoulli function and the azimuthal velocity to vorticity and to the flows.
The rest of this paper is organized as follows: In Sec. II, we briefly recall the derivation of B-H equation, especially the form in spherical coordinate system. In Sec. III, series assumptions for Bernoulli function and azimuthal velocity are listed and briefly discussed. Section IV is about the vorticity decomposition. Sections V-X are devoted to solving different cases of the equation. Section XI offers a summary table and some further discussion.

II. BRAGG-HAWTHORNE EQUATION IN SPHERICAL COORDINATES
For incompressible ideal fluids, the steady motions are described by (steady) Euler equations, which can be written as where v is the velocity and H is the Bernoulli function. LHS (left-hand side) of Eq. (2) is also known as the Lamb vector, denoted as For Eq. (2) to stand, r Â L ¼ 0 is required (and is sufficient, assuming that the flow is in a simply connected domain). A (non-irrotational) flow meeting this requirement is a generalized Beltrami flow 12 (or is further a Beltrami flow if L ¼ 0 in the whole domain 13 ).
In the axisymmetric case, it is possible to define the Stokes stream function w and an independent function C so that v r ¼ 1 r 2 sin h @w @h ; v h ¼ À 1 r sin h @w @r ; in spherical coordinate system ðr; h; uÞ, or in cylindrical coordinate system ðq; /; zÞ.
With the stream function, Eq. (1) is automatically satisfied. Taking the spherical coordinate case as an example, as all derivatives with respect to azimuthal angle u vanish, vorticity is @C @h e r þ À 1 r sin h @C @r e h þ À 1 r sin h @ 2 w @r 2 þ cosh r 3 sin 2 h @w @h À 1 r 3 sin h @ 2 w @h 2 e u (5) and the three components of Lamb vector are L r ¼ À 1 r sin h @w @r À 1 r sin h @ 2 w @r 2 þ cos h r 3 sin 2 h @w @h À 1 r 3 sin h @ 2 w @h 2 À C r sin h À 1 r sin h @C @r ; (6) L h ¼ C r sin h 1 r 2 sin h @C @h À 1 r 2 sin h @w @h Â À 1 r sin h @ 2 w @r 2 þ cos h r 3 sin 2 h @w @h À 1 r 3 sin h @ 2 w @h 2 ; (7) L u ¼ 1 r 3 sin 2 h @w @r @C @h À @w @h @C @r : 8 > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > : As the flow is steady, each of H and C is a function of w only. 14 Thus, we have the following "chain rules" for HðwÞ and CðwÞ: In this case, the Lamb vector components in Eqs. (6)-(8) simplify to L u ¼ 0: Equation (2), if rewritten with Eq. (9) and Eqs. (10)- (12), becomes an identity in e u direction (both sides vanish) and two identical scalar equations in e r and e h direction, and both read This is the Bragg-Hawthorne equation 1 (in spherical coordinate system). It is also referred to as Hicks equation 15 or Squire-Long equation. 16,17 The Grad-Shafranov equation in ideal magnetohydrodynamics also has the same form. 18,19 Changing the variable from h to cosh, we have @w @h ¼ Àsinh @w @ cosh ð Þ and @ 2 w Compared to Eq. (13), LHS of Eq. (14) only has two terms. Considering sin 2 h ¼ 1 À cos 2 h, Eq. (14) has a simple form with respect to the variable cosh. This is convenient. A lot of solutions in this paper are based on this form of the equation.
It is worthwhile to point out that a similar form of Grad-Shafranov equation (using cosh as the variable) is presented by Low and Lou (1990) in a solar magnetic field study paper. 5 The above derivation of B-H equation also applies to the case in cylindrical coordinate system ðq; /; zÞ. In that, Eq. (13) or (14) takes the well-known form as This equation, either in the form of Eqs. (13), (14), or (15), is the governing equation for ideal incompressible axisymmetric steady flows. For a specific case of the functions HðwÞ and CðwÞ, if we can find the stream function w satisfying the equation, velocity can be calculated and Euler equation (2) is solved.
In general, however, this could be difficult. Mathematically, HðwÞ and CðwÞ can be arbitrary (smooth) functions of w. They can make the equation complicated. Nevertheless, when HðwÞ and CðwÞ are in certain forms, the equation can be linear and solvable.

III. LINEAR CASES OF BRAGG-HAWTHORNE EQUATION
LHS of the equation, either in Eqs. (13), (14), or (15), is linear with respect to the stream function w. If the two terms on the RHS are also linear, the equation is linear.
Assuming H 0 , w 0 , k, k 1 , k 2 , a and a 0 , a 1 , a 2 , a 3 are constants, for the first term on RHS to be linear, we have For the second term to be linear, we can have 1 a 3 , a 0 1 2 a 2 1 a 2 ). Mathematically, some constants in the list can be absolved or combined. However, keeping them in the above format will make it convenient for discussion.
Constant w 0 in C 1 , C 2 , and C 4 does not impact the equation and thus has no impact on the solutions. However, when C is interpreted linking to azimuthal velocity, w 0 can bring big differences for azimuthal velocity and for the flow. For C 3 , most of the feasible cases require w 0 ¼ 0. However, to maximize generality, w 0 is kept here as a constant, and the different impacts of w 0 ¼ 0 and w 0 6 ¼ 0 will be discussed specifically.
The cases in the list are not always independent. For example, H 1 can be considered as a special case of H 2 (when k ¼ 0); H 4 is a "combination" of H 2 and H 3 . However, the equations and solutions for these cases are quite different and are worthy to be discussed separately. For this consideration, they are listed as separated cases.
C 3 and C 4 actually give the same form of equation (with different denotations for the constants). So, they would share the same stream function as mathematical solution of the equation. However, when the same solution is respectively interpreted for C 3 and C 4 , the flows can be significantly different (on velocity, vorticity, physical feasibility, etc.). For completeness consideration, they are treated as two separated cases in the list.
In a more general view, for linear cases, when the two terms on RHS are expanded, they can become up to four terms, which can be written on RHS of the following "model equation" as or in cylindrical coordinates K 1 , K 0 , A 1 ; and A 0 here represent certain form of the constants in the list. Different cases of function H and C are associated with different combinations of presentation or absence of these four terms.
Obviously, with H 1 -H 4 and C 1 -C 4 , there are 16 combinations, each making Eqs. (13), (14), or (15) a linear equation of w (i.e., each is related to a combination of the K 1 , K 0 , A 1 ; and A 0 terms in the model equation). For convenience, we will denote these 16 combinations as H m C n (m, n ¼ 1, 2, 3, 4) going forward in this paper.

IV. VORTICITY DECOMPOSITION AND FLOW PROPERTIES
Before exploring solutions of particular cases, it is worthy to further investigate impacts of the two functions, H and C, to vorticity and to the flow.
Recalling Eqs. (3) and (4), we have v u ¼ C rsinh (or v / ¼ C q in cylindrical coordinates). Function C is indeed the linear azimuthal velocity multiplied by distance to the z-axis. 2pC is often considered as circulation along a circle around z-axis. From another angle of view, C can also be considered as the "azimuthal velocity moment" with respect to z-axis. As mass density is a constant in incompressible fluid, C is also representing angular momentum density of the fluid (with respect to z-axis).
In any of these considerations, an assumption (or a restriction) on C is basically an assumption (restriction) on azimuthal velocity related to the stream function w. In other words, C indicates how the fluid is moving around the symmetric axis. It is intuitive to expect that C has close relation with vorticity of the flow.
Bernoulli function H, on the other hand, is the inverse gradient of Lamb vector, which is just the cross-product of velocity and vorticity. It is also natural to expect that H is closely related to vorticity.
With Eqs. (3), (9), and (13), vorticity in Eq. (5) can be rewritten as which is also Same relation exists as well in cylindrical coordinates, where we have Either Eq. (18) or Eq. (19) indicates that vorticity x can be split into two portions. One is in the velocity direction, proportional to velocity with a scale factor dC dw . The second portion is on the azimuthal direction, with magnitude proportional to dH dw (and to the distance to z-axis).
In such a sense, it can be considered that functions H and C construct vorticity, or more specifically, derivatives of H and C compose vorticity in the flow. dC dw decides vorticity in the direction of velocity; dH dw decides vorticity in azimuthal direction. Although these two directions in general are not orthogonal to each other, these two components make up the vorticity at each point in the flow. Equation (18) or Eq. (19) thus can be considered as a decomposition of vorticity. It applies to all flows that B-H equation stands (i.e., all axisymmetric steady flows of incompressible ideal fluid).
Actually, the relations of H and C to the vorticity components are obvious and are often implied in derivation of equations (see, e.g., Batchelor, 14 Sec. 7.5). However, explicitly putting them in the form of Eq. (18) or Eq. (19) straightforwardly shows the impacts of H and C to vorticity and to the flow.

ARTICLE
scitation.org/journal/phf For convenience, we will further denote x B dC dw v, x A Àrsinh dH dw e u (or x A Àq dH dw e / in cylindrical coordinates) and call x B and x A Beltrami vorticity and azimuthal vorticity, respectively. Total vorticity then is sum of these two vectors, i.e., As Beltrami vorticity x B is parallel to velocity, it does not impact Lamb vector. We can just count on the azimuthal vorticity x A when calculating Lamb vector. In other words, regardless of x B , we always This is consistent with the fact that (in ideal incompressible axisymmetric steady flows) Lamb vector is always perpendicular to the azimuthal direction (except for points in the symmetric axis). Following this is also the corollary that Lamb vector is always coplanar with the symmetric axis in such flows.
Applying this decomposition to flows in the list in Sec. III, for the H 1 C n cases, as dH dw ¼ 0, azimuthal vorticity x A vanishes. Total vorticity contains solely the Beltrami part x B . Lamb vector thus vanishes as well. For H 1 C 1 , x B also vanishes, so does the total vorticity. It is hence an irrotational/potential axisymmetric flow. For H 1 C 2 , Vorticity is parallel to velocity with a nonconstant coefficient. This is a non-linear Beltrami flow (i.e., the velocity field is a non-linear Beltrami field). For Vorticity is parallel to velocity with a constant coefficient. This is a linear Beltrami case. H 1 C 4 , on the other hand, is a non-linear Beltrami flow again, as In the H 2 C n family, azimuthal vorticity x A ¼ Àkrsinhe u (or x A ¼ Àkqe / in cylindrical coordinates). It is proportional to the distance to z-axis. For H 2 C 1 , x B vanishes. Total vorticity is made up solely by x A and thus is proportional to the distance to z-axis as well. This is actually the assumption of the Fraenkel-Norbury family of vortex rings 10,11 (including the Hill's vortex without swirl). As can be expected, the general solution of H 2 C 1 will include Hill's vortex as a special case. H 2 C 3 is adding the Beltrami vorticity part x B ¼ av to H 2 C 1 . This is actually the case of Hill's vortex with swirl, which will be discussed in Sec. IX. For H 2 C 2 and H 2 C 4 , both Beltrami vorticity x B and azimuthal vorticity x A present in the total vorticity, and the flows have more complicated dependency to H and C.
In the H 3 C n and H 4 C n cases, as dH dw explicitly contains w, azimuthal vorticity is now directly tied with w (in opposite to H 1 C n or H 2 C n where w is not explicitly involved in the expression of azimuthal vorticity). The equation is still linear, but this azimuthal vorticity brings an extra term with the unknown function w to the equation and makes the solutions complicated. These cases will be further discussed in Sec. X.
As an aside, the second term on RHS of the B-H equation, C dC dw , has a notable property of symmetric. This can also be analyzed by applying the vorticity decomposition.
Mathematically, function C can change to the opposite sign without impacting the equation. If w is a solution of the B-H equation related to CðwÞ, w is also a solution of the equation related to ÀCðwÞ. Physically, this implies that flows governed by B-H equation are "twoway flows." The azimuthal velocity v u ¼ C rsinh or v / ¼ C q can change to the opposite direction and the flow still satisfies the same equation. Equivalently speaking, for each flow as a solution of B-H equation, there exists its chiral symmetric flow (with opposite azimuthal velocity) to form a pair.
From vorticity decomposition point of view, when CðwÞ changes sign (and w remains the same), velocity undergoes a half-round rotation along e u direction without changing chirality [i.e., from changes to a new vector that, in general, is neither rotational symmetric nor chirally symmetric to the original vorticity before C changes sign. Despite the change on vorticity, however, Lamb vector remains the same, and so does the equation (and the solution).

V. AXISYMMETRIC POTENTIAL FLOW (H 1 C 1 )
Combination H 1 C 1 is the simplest one in the list. As discussed in Sec. IV, all flows in this case are axisymmetric potential flow.

As both terms on RHS of B-H equation vanish, Eq. (14) becomes
Let w ¼ RðrÞHðcoshÞ and apply variable separation method. With the separation constant denoted as nðn þ 1Þ, the two separated questions are Solution of Eq. (21) is R ¼ k 1 r nþ1 þ k 2 r Àn . k 1 , k 2 , n here, as well as k 3 , k 4 , k 5 , k 6 , k 0 , k; and K n in expressions later in this paper are independent constants. Mathematically, they can be arbitrary real numbers, or theoretically even be imaginary numbers in some cases. However, to have the stream function and velocity physically meaningful, restrictions may apply to them.
In this paper, we will consider a flow (from a solution of the B-H equation) with no singular velocity (i.e., no infinite nor discontinuous velocity) in a domain as "physically feasible." When the domain is not the whole space, the boundary conditions can be provided by a solid boundary or by other flows outside the domain. In the latter case, the flows inside and outside the domain should have zero normal velocity and continuous tangential velocity at the boundary. Equivalently, this requires inside and outside flows to have constant (normally zero) stream function and continuous stream function gradient on the boundary.
For Eq. (22), H is a function of cosh and derivative € H is with respect to cosh. If we replace the unknown function HðcoshÞ with sinhTðcos hÞ and apply the identities dðsin hÞ

Physics of Fluids
Combining it with solution of Eq. (21) gives the general solution of Eq. (20) as With special settings for the constants, this solution can give specific solutions in simple forms. For example, when n ¼ 1, If k 1 and k 2 are both non-zero, the radial portion in Eq. (24), In the case that such r 0 is a positive real number, stream function w is zero at the surface of the sphere r ¼ r 0 . The flow outside the sphere is a potential flow around that sphere (while the flow inside the sphere is singular at r ¼ 0 and thus is not physically feasible).
In this case, if we further have n ¼ 1, . This is the well-known potential flow past the sphere r ¼ a.
When n > 1, the associated Legendre functions in Eq. (24) can give more complicated potential flows around the sphere. Some examples can be found in Sec. VII [i.e., Eq. (37b) and Fig. 5].
As a property of the C 1 flows (i.e., flows in any H m C 1 case), as long as function C is a constant, its value does not impact the equation and thus will not impact the solution. Azimuthal velocity v u ¼ C r sin h hence is independent of the solution. In principle, velocity ðv r ; v h Þ in the meridian plane derived from solution of an H m C 1 equation can be with any azimuthal velocity v u ¼ C r sin h as long as C is a constant. When C 6 ¼ 0, this v u introduces a circular movement with uniform angular momentum density with respect to z-axis. This is an "irrotational rotation" (as it does not impact the vorticity). In this sense, constant C brings a rotation to H m C 1 flows without impacting vorticity nor impacting velocity in the meridian plane.
Such a rotation has singular velocity at z-axis, though. It is physically feasible only for special cases (e.g., in a domain that is not overlapping z-axis. An example on H 2 C 1 can be found in Sec. VIII, which is the swirled case of Fraenkel-Norbury solutions.).
Case H 1 C 1 can also be solved in cylindrical coordinates. Equation (15) for H 1 C 1 is Let w q; z ð Þ ¼ PðqÞZðzÞ, Eq. (25) can be separated to two equations as q € P À _ P þ kqP¼ 0 and € Z À kZ ¼ 0 (k is the separation constant). The latter one is solved by Z The former one can be converted to a Bessel equation of FðxÞ (of order 1) by replacing PðqÞ by xFðxÞ (with x ¼ ffiffi ffi k p qÞ and thus is solved by linear combination of the Bessel function of the first and second kind, J 1 and Y 1 . Mathematically, the general solution for Eq. (25) is This solution is actually obtained and published previously (1963) by Berker. 20 When is unbounded when z ! þ1 and/or z ! À1, and thus a feasible flow can only exist in a domain excluding upper or lower "end" of the z-axis. When k < 0, with proper k 3 and k 4 , Z z ð Þ is periodic (and bounded) along the whole z-axis. (Alternatively, in this case equation , which is equivalent to the above solution in exponential function form, and is periodic along zaxis.). However, in this case, ffiffi ffi k p is not a real number, and the two Bessel functions in Eq. (26) become modified Bessel functions (of the first and second kind, respectively), which are unbounded on z-axis (i.e., q ¼ 0) or when q approaches infinity.
Overall, physical feasibility of Eq. (26) is limited and needs to be carefully considered according to domain of the flow and boundary conditions, etc.
In passing, w ¼ k 1 q 2 þ k 2 z þ k 3 is also a solution of Eq. (25). This is not from variable separation method and thus is not included in Eq. (26). When k 2 ¼ 0 and k 3 ¼ 0, it converts to w ¼ k 1 q 2 , which is equivalent to the special form of Eq. (24), i.e., the uniform flow w ¼ k 1 r 2 sin 2 h.

VI. AXISYMMETRIC NON-LINEAR BELTRAMI FLOW (H 1 C 2 )
The solution of H 1 C 1 can be "extended" to solve H 1 C 2 . Taking the spherical coordinates case first, Eq. (14) for H 1 C 2 is This can be considered as the homogeneous equation Eq. (20) plus a constant inhomogeneous term on the RHS. It thus can be solved by adding a particular solution, w Ã ¼ À a 4 r 2 to general solution of Eq. (20). This gives Compared to Eq. (24), the particular solution w Ã brings new features for the flow of Eq. (28). It actually turns the potential flow in Eq. (24) into a non-linear Beltrami flow.
If we calculate the velocity from Eq. (28) and compare it with velocity from Eq. (24), w Ã does not impact v r , but it brings an addition v Ã h ¼ a 2 sin h to the polar velocity v h . In the meridian plane, this v Ã h represents a "rotation" around the origin point, with different magnitudes at different polar angles (and is singular on z-axis). Figure 1 shows this v Ã h in x-z plane (when a ¼ 1). The impact of w Ã to azimuthal velocity is more complicated. As is a non-linear function of w, the impact of w Ã to v u is not simply a linear addition.
Theoretically, Eq. (28) can be specified to various flows depending on the constants. However, most (if not all) of them may inherit the singularity from v Ã h and thus need to be within a domain excluding Physics of Fluids ARTICLE scitation.org/journal/phf z-axis. Besides, for v u to be real, it is required that w 0 þ a 2 w ! 0. This may add more restrictions to the flow.
As an example of the feasible cases, when k 1 ¼ a 2 and n ¼ 1, . This is a combination of the uniform flow of Eq. (24) with U ¼ a in z-direction (i.e., w ¼ a 2 r 2 sin 2 h) and the rotation of v Ã h shown in Fig. 1 (that is brought by the particular solution w Ã ¼ À a 4 r 2 ).

Velocity in this case is
q rsinh (as defined in Sec. III, a ¼ a 2 1 a 2 , w 0 is a constant). A real v u requires w 0 À a 2 1 a 2 2 4 r 2 cos 2h ð Þ ! 0. When w 0 < 0, this requirement is satisfied outside of the revolution-solid defined by hyperboloid r 2 cosð2hÞ ¼ 4w 0 . This is the domain the flow exists (as shown in Fig. 2). Velocity has no singularity in this domain. On the surface of the hyperboloid, w is a constant, v u is zero, and velocity is tangent to the surface. The requirement of w 0 < 0 is critical to ensure feasibility of the flow. If w 0 > 0, the hyperboloid has two sheets. The domain [in which w 0 À a 2 1 a 2 2 4 r 2 cos 2h ð Þ ! 0] intersects z-axis. In that z-axis segment, v h and v u are singular and the flow is not feasible.
Similar to the spherical coordinate case, H 1 C 2 in cylindrical coordinates can be solved by adding a particular solution w Ã ¼ À a 4 z 2 or

VII. BELTRAMI SPHERICAL VORTEX (H 1 C 3 )
As discussed in Sec. IV, for H 1 C 3 , as x ¼ av, the velocity field is a linear Beltrami field. The velocity in this case is an eigenvector of the curl operator (and a is the eigenvalue).
In Refs. 21 or 22, a solution of velocity field for this case is given by directly solving x ¼ v with variable separation method (i.e., without employing stream function). As flows in this case have some Physics of Fluids ARTICLE scitation.org/journal/phf significant properties and are important for further discussion, here we will re-solve it in stream function form.
General case of C 3 is C ¼ aw þ w 0 . For convenience, consider w 0 ¼ 0 first. Equation (14) in this case is Let w ¼ RðrÞHðcoshÞ, Eq. (30) separates to the following two equations [with nðn þ 1Þ being the separation constant]: If we set a ¼ 1, k 2 ¼ 0, and k 3 ¼ 1, k 4 ¼ 0, consider the case when n is a positive integer, and re-denote constant k 1 as ÀK n , we have the special case of Eq. (33) as Submitting this into Eq. (3), velocity can be found the same as in Refs. 21 and 22, v r ¼ K n n n þ 1 ð Þr À3=2 J nþ1=2 r ð ÞP n cos h ð Þ; This is a family of multi-layer spherical vortices, indexed by n (and will be referred to as Beltrami spherical vortices in this paper). As descripted in Ref. 22, the field is split by zeros of J nþ1=2 ðrÞ (i.e., zeros of v r and v u ) into homocentric spherical layers. Inside each layer, there are n count of vortex rings, separated by the surfaces P 1 n ðcoshÞ ¼ 0 (i.e., v h ¼ 0, v / ¼ 0). Figure 3 shows contours of the stream function when n is 1, 2, and 6. Examples of velocity field (when n ¼ 1; 2; 3) can also be found in Sec. 15 of Ref. 22.
As shown by the contours, the zero surfaces of J nþ1=2 r ð Þ and P 1 n cosh ð Þ (i.e., zero surfaces of w) split the whole space into axisymmetric and coaxial cells that each contains one vortex ring. On these zero surfaces, velocity is only on the tangential direction. Hence, the fluid inside each cell is contained within that cell all the time.
From Eqs. (34) and (35), stream function and velocity have no singularity in the whole space. Thus, such a Beltrami vortex theoretically can steadily exist by itself in the whole space. In other words, it does not have to be moving in a potential flow (like Hill's vortex) or to be surrounded by other flow outside a restricted domain to be physically feasible.
Nevertheless, it is also possible to assemble a Beltrami vortex with potential flow outside a sphere, similar to the case of Hill's vortex.
For the first order vortices in the family (i.e., n ¼ 1), at the spherical interfaces between layers, where J nþ1=2 r ð Þ ¼ 0; thus w ¼ 0, both v r and v u vanish, and the polar velocity v h is proportional to P 1 1 cosh ð Þ¼ Àsinh. This exactly matches the case of irrotational flow past a sphere. Thus, we can have this Beltrami vortex inside the sphere and match it with that irrotational flow outside.
Specifically, in this case, the stream function can be written as Obviously, in this case, the inside and outside flows both have vanished stream function on the surface of sphere r ¼ D, and gradient of stream function is continuous on that surface.
Depending on which solution of r ¼ tan r parameter D is set to, the vortex inside the sphere can have single or multiple layers. Figure 4 shows two examples of such vortices with 1 layer (D % 4:49Þ and 3 layers (D % 10:90Þ, respectively.

ARTICLE scitation.org/journal/phf
When n > 1, v r and v u still vanish on the interfaces between spherical layers, but the tangential velocity v h at these interfaces has more complicated dependency to h (rather than just being proportional to sinh when n ¼ 1). Even in this case, theoretically it is still possible to match the vortex inside a sphere by an (high-order) irrotational flow outside.
Note that for H 1 C 1 , stream function (24) has the same polar angle section as in Eq. (33). With constants properly selected, a highorder (i.e., n > 1) axisymmetric potential flow of Eq. (24) outside a sphere can match the same order of Beltrami vortex of Eq. (33) inside that sphere with zero stream function and continuous stream function gradient at the interface. This is an extension of the "assembled flow" of Eq. (36).
In this case, the stream function is w ¼ K n r 1=2 J nþ1=2 r ð Þsin hP 1 n cos h ð Þ; ½r D; D is one of the zeros of J nþ1=2 r ð Þ; (37a) When w 0 6 ¼ 0 in C ¼ aw þ w 0 , the equation for H 1 C 3 is @ 2 w @r 2 þ sin 2 h r 2 @ 2 w @ cosh ð Þ 2 ¼ Àa 2 w þ aw 0 . This is Eq. (30) plus a constant inhomogeneous term aw 0 on the RHS. It is solved by adding the particular solution w Ã ¼ w 0 a to Eq. (33). (Another way to solve this equation is to replace w by w À w 0 a so to convert it to the same form as Eq. (30). This gives the same result.) In this case, the particular solution w Ã ¼ À w 0 a brings to the flow an additional azimuthal velocity, which is singular at z-axis. Such an H 1 C 3 flow with w 0 6 ¼ 0 is feasible only when it is inside a domain not overlapping z-axis.
In cylindrical coordinates, when w 0 ¼ 0, Eq. (15) for H 1 C 3 , is This can be solved the same way as that for Eq. (25). The solution is in the same form as Eq. (26), with constant k in the PðqÞ or in ZðzÞ section of Eq. (26) biased by a 2 . That can be written as Same as in the spherical coordinates case, adding w Ã ¼ À w 0 a to Eq.

VIII. HILL'S VORTEX AND EXTENSION (H 2 C 1 )
For Hill's vortex, 7 or a bigger group, the Fraenkel-Norbury family of vortex rings, 10,11,13 the flow in the "core region" is featured by (a) swirl free and (b) vorticity is proportional to distance to the symmetry axis. In terms of B-H equation, feature (a) requires v u ¼ 0, and thus C ¼ 0. Considering the vorticity decomposition, vorticity by Eq. (18) becomes x ¼ Àrsinh dH dw e u (or x ¼ Àq dH dw e / in cylindrical coordinates). In this case, feature (b) translates to dH dw ¼ constant. Flows meeting (a) and (b) hence can be described by H ¼ H 0 þ kw and C ¼ 0. Obviously, this is a subset of H 2 C 1 .
As discussed in Sec. V, for any H m C 1 case, the value of constant C is independent of the equation. Mathematically, a solution in such case is compatible with arbitrary constant C. A non-zero C brings non-zero v u which is an irrotational rotation around the z-axis (and is singular in the z-axis, thus is feasible only in a region not intersecting z-axis). Theoretically speaking (i.e., regardless of the singularity), feature (b) does not have to be paired with feature (a) in a flow.
For any flow of H 2 C 1 , (b) is always valid. In this sense, the inner flow in Fraenkel-Norbury family can be considered as the swirl-free subset of H 2 C 1 (while Hill's vortex is a special case in the Fraenkel-Norbury family).
Equation (14) for H 2 C 1 is This is Eq. (20) plus an inhomogeneous term on the RHS. A particular solution for Eq. (40) is w Ã ¼ k 10 r 4 sin 2 h þ w 0 . Adding it to the solution of Eq. (20) gives the solution of Eq. (40) as Similar to the case of H 1 C 2 , per Eq. (3), velocity v r and v h are linear with respect to the stream function w, thus the particular solution w Ã in Eq. (41) is linearly adding an extra velocity v Ã in meridian plane to the flow of Eq. (24). This extra velocity can be calculated as v Ã r ¼ k 5 r 2 cosh, v Ã h ¼ À 2k 5 r 2 sinh. Figure 6 shows such a v Ã when k ¼ 1.

ARTICLE scitation.org/journal/phf
It is this extra velocity that turns the potential flows of Eq. (24) into H 2 C 1 flows whose vorticity is proportional to distance to z-axis. By itself v Ã is unbounded in far field. As vorticity is proportional to distance to z-axis, vorticity is also unbounded when q ¼ rsinh approaches infinity. To avoid such infinite velocity and vorticity, a physically feasible H 2 C 1 flow in general should be inside a bounded region, with proper boundary conditions provided by another flow outside (or by a solid boundary). This is the case of the Fraenkel-Norbury family of vortex solutions.
To better discuss, we can consider a more general situation that an H 2 C 1 flow [i.e., solution (41)] exists inside an axisymmetric domain, A, that (i) on the boundary @A the stream function w vanishes, (ii) outside A (i.e., in domain R 3 À A) exists an axisymmetric flow for which w and rw are continuous with that of the inside flow at @A, and (iii) the inside flow has no singularity in A and the outside flow has no singularity in R 3 À A (including r ¼ 1).
Obviously, such an "assembled flow" is a candidate for the Fraenkel-Norbury vortex ring solution. In other words, we can consider a flow meeting (i)/(ii)/(iii) as an extended or generalized Fraenkel-Norbury vortex solution. In this case, the inner flow does not have to be swirl-free, and the outside flow does not have to be vorticity-free. These make it more general than the original Fraenkel-Norbury solutions.
In general, such a flow could be complicated. However, when n ¼ 1, Eq. (41) can be in a simple form and we have a chance to study the inside flow explicitly.
To avoid singularity, we can further set k 2 ¼ 0, k 3 ¼ 1; and This is a combination of a uniform flow with U ¼ À2k 1 in zdirection [i.e., w ¼ Àk 1 r 2 sin 2 h as a special case of Eq. (24)] and the extra velocity v Ã as shown in Fig. 6. The polar angle factor, sin 2 h; is the same as in the axisymmetric potential flow Eq. (24) when also with k 4 ¼ 0 and n ¼ 1 (i.e., the potential flow past a sphere).
Boundary of the region, @A, is defined by w ¼ 0. This explicitly gives the equation of @A as r 2 k 1 À k 10 r 2 sin 2 h ¼ w 0 . This definition is equivalent to the definition by Fraenkel and Norbury in Refs. 10 and 11, which appears as w ¼ k (k is a constant). As w 0 is a free constant in the solution, w 0 can be set to zero in Eq. (42) and then w is a constant on @A. In this sense, w 0 here is equivalent to Àk in Refs. 10 or 11. By this equation, the shape of @A is significantly impacted by w 0 . Without losing generality, assume k ¼ 1 and k 1 ¼ 1. Figure 7 shows some cross sections of @A with different values of w 0 .
Obviously, these cross sections are consistent in appearance with the numerical results that Norbury presents in Ref. 11.
When w 0 ¼ 0, @A is the surface of the sphere r ¼ ffiffiffiffiffiffi ffi 10k1 k q . When w 0 > 0, @A becomes a "donut shape" inside the sphere. When 2k , @A reaches the extreme case of a "thin core circle" in x-y plane, as defined by r ¼ The case w 0 < 0 is not feasible as @A in this case (shown by the dotted line in Fig. 7) encloses the whole z-axis and velocity is unbounded in it.
When 0 < w 0 < 5k 2 1 2k , as the core region is totally avoiding z-axis, the flow inside the core region can have any irrotational rotation (i.e., with any non-zero C and azimuthal velocity v u ¼ C rsinh ). This is a swirled case of the (generalized) Fraenkel-Norbury solutions.
The case of w 0 ¼ 0 is Hill's vortex. With w 0 ¼ 0, if we further denote k ¼ A and k 1 ¼ A 10 a 2 , Eq. (42) becomes w ¼ A 10 r 2 sin 2 h r 2 À a 2 ð Þ , which is the well-known stream function of Hill's vortex (in spherical coordinates).
In this case, as a segment of z-axis is involved, there should be no irrotational rotation in the core region (so to avoid singular velocity). In other words, Hill's vortex (as the solution of H 2 C 1 with w 0 ¼ 0) is always swirl-free. In comparison, the case of Hill's vortex with swirl, found by Moffat, 8 is actually a solution of H 2 C 3 , which will be discussed in Sec. IX. In that case, total vorticity is no longer proportional to the distance to z-axis, and thus it does not belong to the Fraenkel-Norbury family of vortex rings.
When w 0 > 0, the outside flow with w and rw matching the inside flow at the donut-shape boundary could be complicated. However, when w 0 ¼ 0, as @A is a sphere surface with polar angle factor sin 2 h for stream function, a potential flow past that sphere matches well the inside flow. That forms the well-known case of Hill's vortex surrounded by a potential flow.
Recall that the first order Beltrami vortex Eq. (36a) also has polar angle factor sin 2 h on the layer interfaces. It is also possible to match Hill's vortex at @A with a "hollow" Beltrami vortex outside. In that case, the stream function can be written as

Physics of Fluids
An example of such assembled vortices is shown in Fig. 8 (with K 1 ¼ 1 and D % 4:49). Inside the interface (i.e., sphere r ¼ D, shown by the dark line in Fig. 8) is Hill's vortex; outside the interface is the "hollow" Beltrami spherical vortex.
The stream function contours appear to be similar to that of the first order Beltrami spherical vortex in Fig. 3. However, certain properties are quite different inside and outside the interface. Vorticity inside is only in azimuthal direction, with magnitude proportional to distance to the z-axis; vorticity outside is parallel to velocity at each point, with magnitude proportional to magnitude of velocity. On the interface, outside flow has vorticity (and velocity) tangent to the sphere surface (i.e., in the azimuthal direction only). This matches vorticity of inside flow and provides a continuous vorticity at the interface. Another obvious difference between the inside and the outside flows is that one is swirl-free and the other one is with swirl.
Switching to cylindrical coordinates, equation for H 2 C 1 , is the inhomogeneous case of Eq. (25). Adding a particular solution, w Ã ¼ k 2 4k5þk6 ð Þ q 2 k 5 q 2 þ k 6 z 2 ð Þ þ w 0 to Eq. (26) gives the solution of Eq. (44) as When w 0 ¼ 0 and k 5 ¼ 1, k 6 ¼ 1, the particular solution becomes w Ã ¼ k 10 q 2 ðq 2 þ z 2 Þ. If we consider the special solution of Eq. (25), w ¼ k 1 q 2 [see the note in Sec. V following Eq. (26)] and add to it to w Ã , we can also get the Hill's vortex solution w ¼ A 10 q 2 a 2 À q 2 À z 2 ð Þ by re-denoting k ¼ ÀA and k 1 ¼ A 10 a 2 . In the case that k 5 and k 6 are both positive and k 5 6 ¼ k 6 , w Ã combined with w ¼ k 1 q 2 gives a specific solution w ¼ A 10 q 2 a 2 À k 5 q 2 ð ). This can be considered as an "ellipsoidal shape version" of Hill's vortex.
For H 2 C 2 , another inhomogeneous term, À a 2 , is added to RHS of the H 2 C 1 equation, i.e., Eq. (40) or Eq. (44). In spherical coordinates that leads to and in cylindrical coordinates, that is, These are still inhomogeneous equation of Eqs. (20) and (25), respectively. The particular solutions for Eq. (46) can be found by combining the particular solution of Eq. (27) and the particular solution of Eq. (40). That gives Similarly, particular solution for Eq. (47) can be found by combining the particular solution of Eq. (29) and the particular solution of Eq. (44). That yields Adding Eqs. (48) and (49), respectively, to Eqs. (24) and (26) gives the general solutions of Eqs. (46) and (47).
For H 2 C 2 , vorticity is neither proportional to distance to z-axis, nor in parallel and proportional to velocity. According to decomposition equations (18) and (19), vorticity in this case has a non-zero Beltrami vorticity part and a non-zero azimuthal vorticity part. Thus, it is in a combined direction of azimuthal direction and direction of the velocity.
Compared to the H 1 C 2 or H 2 C 1 case, the restriction in H 2 C 2 for velocity and vorticity to be bounded and for C ¼ a 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi w 0 þ a 2 w p to be real is tighter. A feasible flow in this case would be in a more restricted region compared to that of the H 1 C 2 or H 2 C 1 flow.

IX. HILL'S VORTEX WITH SWIRL (H 2 C 3 )
The Hill's vortex with swirl, discovered by Moffatt 8 and also by Hicks, 15 is a solution of the B-H equation with C ¼ aw, H ¼ H 0 þ kw. This is the case of H 2 C 3 (with w 0 ¼ 0 in C 3 ).

Physics of Fluids
Obviously, this is the inhomogeneous case of Eq. (30), i.e., the equation of H 1 C 3 . A particular solution of Eq. (50) is w Ã ¼ k a 2 r 2 sin 2 h þ w 0 a . Adding it to Eq. (33) brings the solution of Eq. (50) as The same procedure can be employed to solve the H 2 C 3 equation in cylindrical coordinates, With general solution (39) and the particular solution Alternatively, there is a more straightforward approach to solve Eq. (50). It can be described by the following statement (as a theorem): if w Ã is solution of Eq. (20), andw is solution of Eq. (30), then w ¼w þ w Ã is the solution of the following equation: The proof is straightforward. As w Ã is the solution of Eq. (20), w ¼ w Ã is a particular solution of Eq. (54). Whilew is the solution of Eq. (30), which is the associated homogeneous equation of Eq. (54), w ¼w þ w Ã solves Eq. (54).
This theorem is valid for the cylindrical coordinates case as well. In that case w Ã is the solution of Eq. (25),w is the solution of Eq. (38), and w ¼w þ w Ã solves the corresponding equation of Eq. (54) in cylindrical coordinates, which is Solution of Eqs. (20) and (25) can be in many different forms. Each of them will bring an equation in the form of Eq. (54) or Eq. (55) (not necessarily a B-H equation). Theoretically, this theorem enables an approach to solve these equations. Equation (50) happens to be a specific case of Eq. (54). When we take w Ã ¼ k a 2 r 2 sin 2 h þ w 0 a as a solution of Eq. (20) and apply the theorem, we obtain Eq. (51). Similarly, Eq. (52) is a specific case of Eq. (55) when we take w Ã ¼ k a 2 q 2 þ w 0 a as a solution of Eq. (25). Applying the theorem leads to Eq. (53).
As discussed in Sec. VII, solution (33) represents a family of Beltrami spherical vortices. w Ã in Eq. (51) is adding to them an extra velocity v Ã , which by Eq. (3) can be calculated as In a meridian plane, v Ã r and v Ã h form a constant velocity U ¼ 2k in z direction. On the other hand, azimuthal velocity v Ã u has two terms. The first term k a 2 rsinh is related to a rigid rotation around z-axis; the second term is the irrotational rotation discussed previously, which has singularity on the z-axis. As the flows discussed in this case are mainly vortices centered at r ¼ 0, to avoid singularity in z-axis, we will only consider the case w 0 ¼ 0 (i.e., the case without irrotational rotation for v Ã u ) hereafter in this section. In this case, the extra velocity v Ã represents a special spiral movement: uniform velocity in z-direction and a rigid rotation around zaxis. Adding this to the linear Beltrami flows of Eq. (33) yields the H 2 C 3 flows in Eq. (51).
As nature of rigid rotations, azimuthal velocity v Ã u is unbounded in far field. As a result of that, a feasible H 2 C 3 flow should be inside a bounded region with proper boundary conditions. When n ¼ 1 (and k 4 ¼ 0), Eq. (51) has a polar angle factor sin 2 h. This is the same as the irrotational flow past a sphere. Thus, that irrotational flow can provide the boundary conditions. Taking Eq. (57) as the special case of Eq. (51) and matching it with the irrotational flow outside, with some constants re-denoted, the stream function can be written as Note that sinðarÞ ar À cosðarÞ þ k Aa 2 r 2 ¼ 0 can have multiple positive solutions (depending on value of k Aa 2 ). Each solution represents a spherical interface of a closed layer. Figure 9(a) shows the flow of Eq. (58a) in the whole space when A ¼ 1, a ¼ 1; and k ¼ 0:01. Two closed layers exist near the center in this case. Outside of the second layer, the flow is no longer contained in close layers, and velocity increases (unboundedly) as r increases.
The flow inside any closed interface can be matched by a potential flow outside. In other words, it does not have to be the first interface from the center. In the case that D in Eq. (58a) is set to a solution other than the first one from the center, the inner vortex has multiple layers. Figure 9(b) shows the central section of Fig. 9(a) with two layers (D % 7:18) surrounded by potential flow of Eq. (58b). Similar to Hill's vortex without swirl in Eq. (43), the inner vortex with a polar angle factor sin 2 h can also be matched by a hollow Beltrami Spherical vertex (similar to that in Fig. 8). In this case, the inner flow is still the same as Eq. (58a), and the outer flow is from Eq. (33) with dedicated constants to have stream function and its gradient matching the inner flow at the interface, An example of such vortices is shown in Fig. 10 (with A ¼ 1, a ¼ 1; and k ¼ 0:01, D % 7:179 8, b % 1:076 0). Similar as for H 1 C 4 /H 1 C 3 , the H 2 C 4 case has the same equation as H 2 C 3 and thus shares with H 2 C 3 the same solution (51) [or solution (53) in cylindrical coordinates]. Feasibility and properties of the H 2 C 4 flows, however, are quite different from that of H 2 C 3 , even with the same stream function.

X. OTHER AXISYMMETRIC FLOWS (H 3 C N AND H 4 C N )
Some physical background of the H 3 C n /H 4 C n cases can be found in, e.g., Refs. 3, 6, and 23.
For these two families, as H is a quadratic function of w, the first term on RHS, r 2 sin 2 h dH dw or q 2 dH dw , explicitly contains w. Moreover, in the spherical coordinates case, two coordinate variables, r and h, are explicitly coupled with w in this term. This makes it difficult for variables to be separated. In cylindrical coordinates, however, this term involves only one coordinate variable q. This gives a chance for variable separation method to apply.
In the rest of this section, we will concentrate on possible solutions for these two families in cylindrical coordinates. The approach to solve them in spherical coordinates is yet to be explored.

Physics of Fluids
ARTICLE scitation.org/journal/phf to be a combination of Beltrami vortex and a special rotation, and can be extended to multi-layer cases. The other way the new general solutions in this paper can contribute is to provide new flows as explicit analytical solutions of B-H equation. For this, the non-linear Beltrami flow discussed in Sec. VI is an example. For convenience, a summary is given in Table I for different cases and the solutions, as well as typical flow properties, when available.
In the table, presentation and absence of the four terms on RHS of the model equation are indicated by "1" and "0" on the "RHS terms" column. w Ã is representing the particular solution for each inhomogeneous equation. As w 0 in C 3 has big impacts on the equation and solution for H 1 C 3 , cases of w 0 ¼ 0 and w 0 6 ¼ 0 are listed separately. Same are the H 1 C 4 cases with a 0 ¼ 0 and a 0 6 ¼ 0. By nature of the variable separation method, solutions can be different when variables are separated in different ways (e.g., in different coordinate systems). Thus, other ways to separate variables may find other solutions. The solutions obtained so far in this paper are not expected to cover all possible solutions.
As the major attention in this paper is on general mathematical solutions, only a few obvious specific flows were discussed as examples. By considering all possible values of the constants, there could be many more specific flows to explore. Physical feasibility (and stability) of them, however, is to be carefully investigated, especially for those "assembled vortices."

DATA AVAILABILITY
The data that support the findings of this study are available within the article.