A novel fast simulation technique for axisymmetric PWFA configurations in the blowout regime

In the blowout regime of plasma wakefield acceleration (PWFA), which is the most relevant configuration for current and future applications and experiments, the plasma flow that is excited by the ultra-relativistic drive beam is highly nonlinear. Thus, fast and accurate simulations codes are indispensable tools in the study of this extremely important problem. We have developed a novel algorithm that deals with the propagation of axisymmetric bunches of otherwise arbitrary profile through a cold plasma of uniform density. In contrast to the existing PWFA simulation tools, our code PLEBS (PLasma-Electron Beam Simulations) uses a new computational scheme which ensures that the transverse and longitudinal directions are completely decoupled---a feature which significantly enhances the speed and robustness of the new method. Our numerical results are benchmarked against the QuickPic code and excellent agreement is established between the two approaches. Moreover, our new technique provides a very convenient framework for studying issues such as beam loading and short-range wakefields within the plasma cavity.


I. INTRODUCTION
The technique of plasma wakefield acceleration (PWFA), in which an intense, ultra-relativistic drive beam excites strong accelerating fields as it moves through a dense plasma column, is one of the most promising schemes for achieving the unprecedented acceleration gradients necessary for linear collider or compact free-electron laser (FEL) applications [1,2]. Unlike the early research efforts in this area [3][4][5], most contemporary iterations of this concept are based on the so-called blowout regime [6,7], in which the density of the driver is comparable to (or considerably higher than) that of the plasma background. In this case, the drive beam expels plasma electrons from its path in such a way that a co-moving cavity (or bubble) is created in its wake. The corresponding plasma flow is highly nonlinear, which makes the development of a rigorous, analytical theory for this regime extremely difficult. Most efforts along these lines have focused on a phenomenological treatment of the problem, using insight gained from simulations [8][9][10]. Recent attempts to construct an analytical description of the blowout regime from first principles have been surprisingly successful, but the resulting expressions are only applicable in the limit of small driver charge and small driver dimensions [11]. As a result, simulation codes remain the main tool for treating the PWFA problem in its general form [12][13][14][15].
Moreover, the beam loading effect associated with the presence of a realistic witness bunch and beam instability issues related to short-range wakefields (longitudinal and transverse) induced within the plasma bubble are topics of considerable importance when assessing the feasibility of the PWFA concept [16]. In this paper, we present the outline of a novel PWFA simulation code that, apart from treating the basic problem of electron acceleration, also offers a natural and convenient framework for dealing with more advanced subjects like the ones mentioned above. Starting entirely from first principles, our semi-analytical formalism is based on the assumption that we are only interested in the steady-state regime of the interaction. This quasi-static approximation implies that we neglect the longitudinal plasma non-uniformity and other fast effects associated with the injection process. Another crucial assumption is that the PWFA configuration under consideration is axially symmetric. Restricting our analysis to the 2D case enables us to simplify our treatment while retaining most of the underlying physics. At the same time, we can readily take into account the finite size of both the driver and the witness bunches, which can have an arbitrary density profile.
Our development is organized as follows: In Section II, we formulate the equations for the steadystate electromagnetic field excited by a beam moving with the speed of light in a uniform plasma. Section III then derives the single particle equations of motion for the plasma electrons, while Section IV introduces some basic concepts regarding the description of the plasma flow in terms of macroparticles. The remaining details of the computational algorithm are given in Section V, while the results of a simple numerical study of the beam loading effect are presented in Section VI. As part of the latter, we also include a comparison between our code and an existing particle-in-cell simulation tool. Our analytical derivation is completed in Section VII, where we show how our code can be adapted for the calculation of short-range wakefields inside the plasma cavity. Finally, Section VIII summarizes the main results of this paper.

II. EQUATIONS FOR THE ELECTROMAGNETIC FIELD
In this section, we formulate the equations that describe the plasma dynamics behind a finite-size driver moving through a cold plasma with the speed of light. To start with, we assume that the driver propagates along the z axis in the positive direction. As we have already mentioned, our analysis is restricted to the steady-state case, so sufficient time is assumed to have elapsed since the injection stage and the equilibrium plasma density n 0 is taken as uniform in space (moreover, we disregard the motion of the ions). In order to incorporate the effect of beam loading, we also assume that a finite-size, ultra-relativistic witness beam follows the driver through the plasma column at a fixed distance behind it. Both the driver and the witness bunches are assumed to be on-axis, radially symmetric and non-evolving. Following the standard convention, we normalize the time t to ω −1 p , length to k −1 p , and velocities to the speed of light c. Here, ω p = ck p is the plasma oscillation frequency, given by ω p = 4πn 0 e 2 /m = c √ 4πn 0 r e (r e = e 2 /mc 2 is the classical electron radius). We also normalize momenta to mc, fields (electric and magnetic) to mcω p /e, potentials (scalar and vector) to mc 2 /e, charge densities to n 0 e and current densities to en 0 c.
Assuming that all potentials and fields depend on z and t solely through the combination ξ = ct−z, which expresses the longitudinal position with respect to a fixed point of the driver, we find the following equations for the non-zero components of the (scaled) electric and magnetic fields E r , E z and B θ in the cylindrical coordinate system (for the current density, we use the relation j = −nu, where the negative electron charge is explicitly taken into account): where n is the scaled plasma electron density, n ext = n d + n w is the total "external" density due to the driver and witness beams while v r and v z are the radial and longitudinal components of the electron collective velocity. We also need the equation for the divergence ∇ · E, which yields From Eqs. (1b) and (1c), we can introduce the pseudo-potential ψ such that The function ψ is the difference between the scalar potential φ and the longitudinal component of the vector potential A z , i.e. ψ = φ − A z . Combining Eq. (1a) with (2) , we obtain which gives us the following equation for ψ: Additionally, we have an equation for the gradient ∂ ξ ψ, namely We also need an equation for the azimuthal magnetic field B θ . To start with, we differentiate Eq. (1a) with respect to r and substitute ∂E z /∂r with the value extracted from Eq. (1c). This yields the relation ∂ ∂r We will discuss how Eqs. (5) and (7) are solved numerically in Section V.

III. EQUATIONS OF MOTION FOR THE PLASMA ELECTRONS
In addition to the equations for the electromagnetic field derived in the previous section, we need the equations of motion for the plasma electrons. Using the relation dξ = (1 − v z )dt for the differentials dt and dξ, these equations can be written as where p r = γv r and p z = γv r are the radial and longitudinal components of the momentum vector (note that these are now the individual particle momenta). We also have a corresponding equation for the relativistic γ factor: There is an important integral of motion in this problem, given by γ − p z − ψ = const. This can be shown using Eqs. (8b), (9) and (3) [17]. Assuming that the plasma electrons are at rest at ξ = ξ init , where ξ init is the initial (lower) ξ-value corresponding to the front of the driver, we obtain γ − p z − ψ = 1 or where we have also used the fact that ψ(r, ξ = ξ init ) = 0, a property which will be verified later on. From this, we can express p z and γ in terms of p r and ψ via the relations Next, we re-write the equation for the radial momentum p r as where we have used Eqs. (3) and (10). In addition to the above, we also require an equation for the radial orbit. Starting from dr/dt = v r , we obtain We should point out that the usefulness of ξ = t − z as an independent variable instead of the time t (in the particle equations) is enhanced by the fact that the former is a monotonically increasing function of the latter (recall that dξ/dt = 1 − v z > 0). The single-particle equations of motion should be supplemented by the continuity equation for the plasma density n, namely IV. MACROPARTICLES In our computational algorithm, the plasma electrons are represented as a sum over macroparticles. Each macroparticle is characterized by the dimensionless charge q i , the coordinate r i (ξ) and momenta p zi (ξ) and p ri (ξ). The coordinates and momenta satisfy the equations of motion (12) and (13), with the fields on the right hand side of the equations taken at the positions of the particles. The plasma electron density n and the corresponding currents are represented as a sum over macroparticles according to the relations where v zi and v ri are the macroparticle velocities and the summation goes over all macroparticles in the system. The particular form of the weights in front of the delta functions in Eqs. (15) is chosen in such a way that the continuity equation (14) is automatically satisfied by the expressions given above.
The weights q i are determined by the initial coordinates of the macroparticles, r i (ξ init ). Integrating the first of Eqs. (15) from r 1 to r 2 with the weight r and using v zi (ξ init ) = 0 and n(r, ξ init ) = 1 (both stemming from the fact that the plasma is at rest in front of the driver) we obtain To approximately satisfy this equation we choose q i = δr i r i (ξ init ) where δr i is the separation between the particles; then the sum on the right-hand side becomes i δr i r i (ξ init ) ≈ 1 2 (r 2 2 − r 2 1 ). Of course, to minimize the error stemming from this approximation one has to use many particles with the distance between them δr i much smaller than the transverse size of the plasma flow.

V. CALCULATION OF THE FIELDS
In this section, we focus on the solution of the equations for the plasma electromagnetic fields (5) and (7). To start with, we decompose the magnetic field into two parts via This quantity, which expresses the contribution of the driver and witness beams to the total magnetic field B θ , is given in an analytical fashion by The equation for the part of the magnetic field due to the plasma electron flow then becomes We need to solve this equation with the boundary conditions: B θ = 0 at r = 0 and B θ → 0 when r → ∞.
For the radial current nv r , we use the last of Eqs. (15), which can be written as where , ξ) and we have also made use of Eq. (10). In order to simplify the notation, we have temporarily suppressed the ξ-dependence of the various macroparticle quantities.
For the longitudinal plasma current nv z , we have the relation where we have again used Eq. (10). The radial derivative ∂(nv z )/∂r is easily calculated to be where the prime denotes differentiation with respect to r. Considerably more work is required for obtaining the longitudinal derivative ∂ ξ (nv r ). The end result of the corresponding calculation is where use has been made of Eqs. (12) and (13). This transformation of the derivative ∂(nv r )/∂ξ can be also found in Ref. [18]. Substituting Eqs. (22) and (23) into Eq. (19) forB θ , we observe that on the right-hand side there is a sum of terms each of which is either proportional to δ(r − r i ) or δ (r − r i ). For r i < r < r i+1 (where we have tacitly assumed that the macroparticles have been labeled according to increasing radius), the plasma currents are zero andB θ satisfies the relation ∂ r [(1/r)∂ r (rB θ )] = 0. Thus, the magnetic field between the delta functions with labels i and i + 1 can be represented as where a i and b i are functions of ξ. Our objective now becomes to find the corresponding matching conditions at the locations of the delta functions, r = r i (ξ). These conditions are derived in Appendix A and are formulated as a system of linear equations for a i , b i : where the coefficients A i , B i and C i are Eqs. (25) can also be written in matrix form, Applying the recursion rule of Eq. (27) for 1 ≤ i ≤ N (where N is the number of simulation particles), we can obtain the magnetic field pattern at a given ξ given the pseudo-potential ψ and the macroparticle positions and momenta. Apart from the condition b 0 = 0 (which avoids a singularity at r = 0 and guarantees that B θ = 0 at r = 0), we typically assume that a M = 0, which ensures a 1/r field decay after the last particle, and hence B θ → 0 when r → ∞.
To complete the electrodynamics-related part of our derivation, we now consider Eq. (5). The boundary conditions for this equation are: ψ is finite at r = 0 and ψ → 0 when r → ∞. Substituting the second of Eqs. (15) to the right-hand side of Eq. (5), we see that it contains terms that are proportional to δ(r − r i ). For the region between two neighboring macroparticles (r i < r < r i+1 ), we have zero plasma current so (1/r)∂ r (r∂ r ψ) = −1. This relation yields a vacuum solution of the form ψ =ā i +b i ln r − r 2 /4, whereā i ,b i are functions of ξ. To find relations between the adjacent coefficientsā i ,b i , we integrate Eq. (5) over the radius r from r − i = r i − 0 to r + i = r i + 0 through the delta function δ(r − r i ). We obtain the following matching conditions for the pseudo-potential ψ, With these matching conditions, we find thatā i +b i ln r i =ā i−1 +b i−1 ln r i andb i −b i−1 = q i . For 0 < r < r 1 (i.e. for radii smaller than that of the first macroparticle), we have ψ =ā 0 − r 2 /4 sinceb 0 must be equal to zero in order to avoid a singularity at r = 0. Using the above, we can add potential terms over successive delta functions in order to construct the solution for ψ. The end result is ψ(r, ξ) = ψ(0, ξ) + ∆ψ(r, ξ), where and Here, we have adopted the condition that the pseudo-potential ψ is always zero at the location of the last macroparticle, i.e. ψ(r N (ξ), ξ) = 0. For large enough radial coordinate of the last electron in the system, r N (ξ), this condition approximates the boundary condition ψ → 0 when r → ∞. From Eqs. (29) and (30), one can easily derive the partial derivatives of ψ that are needed for evaluating the A i , B i , C i coefficients discussed earlier. In particular, we have for the radial derivative and ∂ ξ ψ(r, ξ) = ∂ ξ ∆ψ(r, ξ) − ∂ ξ ∆ψ(r N (ξ), ξ) for its longitudinal counterpart, where and we have assumed that dr N (ξ)/dξ = 0. Incidentally, it is worth mentioning that one of the ξderivatives mentioned above is directly related to the calculation of the on-axis longitudinal electric field, one of the main figures of merit for the PWFA. Specifically, the force per unit charge for an onaxis electron is given by F z /e = −E z (0, ξ) = ∂ ξ ψ(0, ξ) = −∂ ξ ∆ψ(r N , ξ), a result which is particularly useful for numerical calculations. Lastly, we point out that combining Eqs. (16) and (31) yields ∂ r ψ(r, ξ init ) = 0, which means that ψ(r, ξ init ) = ψ(r N , ξ init ) = 0. This verifies our earlier claim that the pseudo-potential is zero for all r at the front of the driver. Eqs. (18), (24), (27) and (26), along with Eqs. (29)-(32), form the basic results of our field analysis. Their main feature is that a knowledge of the macroparticle positions and momenta at a given ξ suffices for a complete description of the electromagnetic fields at that particular (relative) longitudinal position. This is also reflected in the absence of any derivatives with respect to ξ on the LHS of Eqs. (19) and (5), a property which prevents the coupling of different ξ-values. As a result, one can readily make use of the single-particle equations of motion (i.e. Eqs. (11)-(13)) in order to propagate the particle positions and momenta to ξ + δξ. To avoid any confusion, we rewrite the latter in our updated (macroparticle) notation as and In our algorithm, we use a fourth-order Runge-Kutta technique to numerically integrate the equations of motion for the plasma electrons. Another important point is that, when calculating discontinuous quantities likeB θ , ∂ r ψ or ∂ ξ ψ at the location of a macroparticle (i.e. at the actual location of the discontinuity), we use the arithmetic mean of the two limiting values. This is simply due to the fact that the simulation macroparticles really represent charged sheets formed by the nonlinear plasma flow. In any event, knowing the new positions and momenta enables us to calculate the fields at the updated location. Thus, a full solution of the steady-state PWFA problem (including details such as the shape of the cavity) can be constructed through this systematic, step-by-step process.
We have implemented our algorithm in the computer code PLEBS (PLasma-Electron Beam Simulations) using the MATLAB programming environment. In the next section, we will compare the results obtained with this code against QuickPic simulations.

VI. BEAM LOADING STUDY
So far, we have kept the density profiles of the driver and witness beams entirely general. Here, we choose to specialize to the case where both beams have a Gaussian profile in the transverse and longitudinal directions. In doing so, we find it most convenient to start from the original quantities and then switch to their dimensionless counterparts. To begin with, we assume that the density of the drive beam is given by where n d0 is the peak volume density, σ rd and σ ξd are (respectively) the transverse and longitudinal rms beam sizes and ξ d is the location of the drive beam centroid. Regarding the origin of ξ, we select ξ init = 0 and require this initial value to correspond to the front (or head) of the drive beam. This, in turn, is virtually guaranteed if we choose ξ d = N 0 σ ξd , where N 0 ≥ 3. The current profile of the drive beam is where I d0 = 2πcen d0 σ 2 rd is the peak current. The total charge contained in the drive beam is Q d = eN d = √ 2πI d0 (σ ξd /c), where N d = n d0 (2π) 3/2 σ 2 rd σ ξd is the total number of electrons. As far as the witness beam is concerned, we have an analogous Gaussian expression, namely The corresponding current profile is given by with I w0 = 2πcen w0 σ 2 rw , Q w = eN w = √ 2πI w0 (σ ξw /c) etc. At this point, we switch back to the dimensionless variables employed in most of this text (we recall that volume densities are scaled with respect to the plasma density n 0 , while lengths are normalized with respect to the plasma skin depth k −1 p ). The scaled drive beam density is where the tildes (temporarily) denote the scaled notation andñ d0 is given byñ d0 = √ 2/π ν d /(σ 2 rdσ ξd ). Here, ν d ≡ N d r e k p is the dimensionless charge of the driver. This important quantity can be shown to be proportional to the ratio of the number of particles in the beam to the number of plasma particles within a sphere of radius k −1 p [11]. Needless to say, an entirely analogous expression can be obtained for the witness beam.
Dropping the tildes from now on, we use the above results and Eq. (18) in order to determine the external field B (0) θ for our case. The resulting analytical expression is With the above relation, we have everything that we need in order to study the beam loading effect for a Gaussian drive/witness beam configuration. For the drive beam, the parameters we have used are σ rd = 2 µm, σ ξd = 20 µm, I d0 = 2 kA (so that Q d = 334 pC) and ξ d = 3σ rd . As far as the witness beam is concerned, we assume that σ rw = 2 µm, σ ξw = 6 µm, ξ w = 130 µm and I w0 = 2.125/2.55 kA (we have examined two different cases, corresponding to a charge of 106/127 pC). The plasma density n 0 is taken to be 7 × 10 16 cm −3 , which leads to a skin depth of k −1 p = 20 µm. For these parameters, the dimensionless charge of the drive beam is ν d = 0.293, while that of the witness bunch is ν w = 0.09/0.11. Moreover, we have n d0 ≈ 23.7 and n w0 ≈ 25.2/30.2. The fact that these scaled density values are much larger than unity is a typical characteristic of the blowout regime of the PWFA.
The current profiles for the drive and witness bunches are illustrated in the right graph of Fig. 1 (for a witness current of 2.125 kA). In the left graph of the same figure, we plot a large number of electron trajectories whose initial radii r i (ξ = 0) are uniformly spaced on the r-axis. More specifically, our simulation run makes use of 2000 macroparticles with 0 < r i (ξ = 0) < 5, though we only plot 400 representative trajectories. The formation of a bubble is evident, with a maximum radius r cav ≈ 1 (about 20 µm) and a length ξ cav ≈ 8 (or 160 µm). In Fig. 2, we plot the on-axis, longitudinal electric field E z (0, ξ) as a function of ξ for the two cases of witness beam current. The solid lines correspond to results from PLEBS, and the dashed lines represent data from QuickPic [19], a popular particle-in-cell PWFA simulation tool. The left graph shows the comparison for the full range of ξ (up to the bubble collapse), while the right graph zooms in on the location of the witness bunch. According to a well-known pattern, the electric field is initially decelerating but changes sign as one moves further away from the driver. Moreover, the beam loading effect only influences the field in the neighborhood of the witness beam. Overall, very good agreement is observed between the two approaches, while our technique appears to offer an advantage in terms of computation speed. In particular, our MATLAB-based code only takes a few minutes to run on a simple desktop environment, whereas QuickPic requires a parallel setup with 128 cores in order to achieve the same performance.

VII. WAKEFIELD CALCULATION
In this section, our aim is to demonstrate that the algorithm which was described in the previous parts of this paper can also be used to calculate the short-range wakefields induced within the plasma bubble. For the basic definitions, notation and analytical results regarding the topic of wakefields in a PWFA context, we rely on work presented elsewhere [20]. To start with, we turn to the problem of the wakefields excited by an on-axis point charge q located inside an axisymmetric plasma cavity at the longitudinal position ξ = ξ 0 . Our goal is to calculate the longitudinal wakefield immediately behind the charge, at ξ = ξ + 0 . We denote byẼ r (r, ξ),Ẽ z (r, ξ) andB θ (r, ξ) the components of the modified field in the presence of the point charge q. As is shown in [20], the transverse field components can be expressed asẼ where E r , E z and B θ are the original components, h(ξ) is the step function equal to one for positive arguments and zero otherwise, ∆E r and ∆B θ denote the change of the field due to the charge q, and the terms with the delta function represent a shock wave characterized by the radial profile D(r, ξ 0 ). The latter quantity satisfies the relation The longitudinal electric field E z also exhibits a discontinuity, which is expressed bỹ The jump ∆E z (r, ξ) represents the longitudinal wake generated by the charge q immediately behind it. It can be related to the D-function by means of the relation Thus, calculating the longitudinal wakefield involves solving Eq. (42) for a given (that is, known from a previous calculation) plasma density n(r, ξ). Using the first of Eqs. (15), we can re-write (42) as where A i has been defined in Eq. (26) and we have used Eq. (10). The above equation can be solved using an iterative technique analogous to the one we employed in solving Eq. (19) (we note the similarity between these two equations, which have essentially identical left-hand sides). Specifically, since the RHS of (45) is zero in the intervals r i < r < r i+1 (for i = 1, . . . , N − 1), we seek solutions of the form D =â i r +b i /r, whereâ i ,b i only depend on ξ. Furthermore, we have D =â 0 r + 2ν q /r for r < r 1 and D =b N /r for r > r N . Here, ν q = qr e k p /e is the scaled charge and we have also taken into account the asymptotic behavior of the shock-like field induced by a point charge in a uniform plasma (E r ∝ ν q /r for r → 0, according to [21]). Lastly, our choice ofâ N = 0 ensures that D → 0 at r → ∞. At the location of each delta function, we have the matching conditions of continuity and derivative jump for D, given by These conditions yield the relationŝ withb 0 = 2ν q andâ N = 0. In matrix notation, the solution of this system is where the transfer matrixM i is given bŷ We note that this particular matrix is also present in Eq. (27). Combining these manipulations, we obtain whereM =M NMN−1 . . .M 1 is a cumulative matrix. Recalling thatb 0 = 2ν q andâ N = 0, we obtain a 0 = −2ν qM12 /M 11 . Combining Eq. (44) with the small-radius expression for D (D =â 0 r + 2ν q /r), we find that the on-axis value of the longitudinal wakefield is given by For the transverse wakefield calculation, the configuration of the leading point charge is somewhat different. In particular, the point charge q is now off-axis, moving with a transverse offset a inside the plasma cavity (we assume |a| k −1 p ). This time, the perturbed electromagnetic field, which now lacks axial symmetry, is given bỹ Assuming that a = ax, the components of the shock profile D (in cylindrical coordinates) are where the radial profile u(r, ξ) satisfies the relation ξ(c/ω p )

VIII. CONCLUSIONS
In this paper, we have developed a novel, steady-state simulation technique that can deal with an axisymmetric PWFA configuration in the blowout regime. In particular, we have studied the propagation of two ultra-relativistic, non-evolving and axially symmetric bunches of arbitrary density profile through a cold plasma of uniform density. After formulating and analyzing the single-particle equations of motion for the plasma electrons, we show how the nonlinear plasma flow of the PWFA can be modeled using simulation macroparticles. In order to obtain a self-consistent description of the interaction, we combine our analysis of the plasma dynamics with the equations that govern the formation of the electromagnetic fields. A crucial feature of our semi-analytical treatment is that, given the macroparticle coordinates at a given longitudinal position within the plasma bubble, we can determine the full pattern of the fields at that particular location. This decoupling between the transverse and longitudinal directions is what makes our approach different from the existing numerical algorithms [13,17,19] and allows us to determine the structure of the bubble in a systematic and efficient way. Using our technique, we have studied the beam loading effect for a Gaussian drive/witness beam configuration. Our results regarding the dimensions of the plasma cavity and the electric field pattern that is established within it are in agreement with those obtained from Quick-Pic, a well-known particle-in-cell PWFA code. In addition to this, our approach appears to offer a relative advantage in terms of computation speed, simplicity and versatility. The latter feature is particularly emphasized by the fact our algorithm can be used directly for the calculation of the short-range wakefields inside the plasma cavity.

IX. ACKNOWLEDGEMENTS
We would like to thank V. Khudik for his important contributions to this work. We are thankful to X. Xu for providing the QuickPic data for the beam loading study and to W. An for the comparison of our algorithm with that of QuickPic. This work was supported by the Department of Energy, contract DE-AC03-76SF00515.
Moving the term c 0 δ(x − x 0 )/x = c 0 δ(x − x 0 )/x 0 to the RHS, we have Since z is still continuous (no delta function derivatives among the driving terms), the added term proportional to z does not affect the matching conditions at x = x 0 . The same is true of the term proportional to the combination z + c 0 h(x − x 0 ), which has a finite discontinuity at x = x 0 . Thus, the sole effect of the change in the form of the LHS of (A1) is the replacement of b with b − c 0 /x 0 in (A5). By reviewing Eqs. (19), (22) and (23), we can easily verify that the equation forB θ is entirely analogous to the one studied above. In particular, the basic analogies are x → r, x 0 → r i and y →B θ , so that we can use Eq. (A5) with and a → A i , b → B i and c 0 → C i . This leads directly to Eq. (25) in the main text.