BBNBI beam divergence and verification

[miekkasarki]

Right now BBNBI ignores the halo component and assumes isotropic divergence using only the horizontal value. All relevant parameters are in the input already, so implementing this is just a few lines of code. (I label this as a bug right now.)

For the verification test, we could have a single beamlet and verify from the beam deposition that the divergence and halo were properly implemented (the distribution should be bivariate double Gaussian).

We could also inject the beam into an uniform plasma and using the mean-free-path for the beam ionization we would have analytical estimate for the beam deposition.

[miekkasarki]

Based on the model in the old BBNBI (in ASCOT4), the beam divergence and halo fraction will be implemented as

$$P(\omega_x,\omega_y)/P_\mathrm{tot}=\frac{f}{\pi X \omega'_h\omega'_v}\left(e^{-\frac{1}{2}(\omega_x/\omega'_v)^2}+e^{-\frac{1}{2}(\omega_y/\omega'_h)^2}\right) + \frac{1-f}{\pi X \omega_h\omega_v}\left(e^{-\frac{1}{2}(\omega_x/\omega_v)^2}+e^{-\frac{1}{2}(\omega_y/\omega_h)^2}\right)$$

Here $f$ is the halo fraction, $\omega_h$ ($\omega'_h$) is the horizontal (halo) divergence, and $\omega_v$ ($\omega'_v$) is the vertical (halo) divergence. This is essentially the same as Eq. (1) here but with the horizontal and vertical divergences separated.

[pietrovincenzi]

Speaking with our DTT NBI design team colleagues, they define halo and core beamlet fractions as in the picture below (numbers are for DTT, obsolete now):

We checked, and this is in agreement with BBNBI definition.
It is important to leave this user-defined parameters:

• halo and core divergences (separately)
• horizontal and vertical divergences
• f, i.e. the power fraction carried by halo

So in total (up to) 5 parameters. For divergences, it would be easy to use the angles corresponding to half of the cones (as it should already be), e.g.:

[miekkasarki]

Great! We already have all those 5 parameters in the input so this is just a matter of implementing that in the code. Shouldn't be an issue. I will write an algorithm that:

1. Draws a random number between [0,1]. If it is smaller than $f$ (the halo fraction) it will be part of the halo.
2. Depending on whether the marker is in halo or not, choose the correct values for the divergences.
3. Samples the Gaussian distributions to get values for $\omega_x$ and $\omega_y$. This is easy since these variables are not correlated.

[miekkasarki]

@rui-coelho I think I found the reason why the beam bends. The momentum in cylindrical basis was not updated. My bad but luckily this issue is not present in the 5.4 release. And it won't be in 5.5.

[rui-coelho]

Is it possible to have it fixed in the main ?

[rui-coelho]

And regarding the definition of the beam divergence, what the Japanese provide as half-divergence angle is the same as in Pietro's figure labeled as "divergence angle" (misleading...). In ASCOT4/old-BBNBI formulae, i'm not so sure if the "divergence" is rather the half divergence angle multiplied by sqrt(2). Dr. Google sent me to some sites on optics and lasers and in all of them it is widely accepted that :

tan(half-divergence angle) = (difference in radii of two beam cross sections along propagating path) / (distance of the two cross sections along propagating path)

And the "radius" of the beam is usually defined as the width at 1/e^2 of the beam intensity. In your definition of the "probability"/power the "radius" of the beam will get 1/e^2 if we use and angle of half-divergence * sqrt(2).

It is thus very important to clarify if:

1. For bbnbi5 what we put in the Injector class object is the half-divergence angle or the full one
2. For bbnbi5 that "angle" implemented results in a beam "radius" that increases with tan(angle) or tan(angle*sqrt(2))

[miekkasarki]

I pushed a fix for the bending issue and implemented the divergence. Also BBNBI is now part of the automatic testing which checks that both the beam divergence is correct and the ionization rate based on how the ionized markers are distributed.

The code needs some polishing but you can already try it out in branch https://github.com/ascot4fusion/ascot5/tree/hotfix/48-bbnbi-beam-divergence-and-verification

The divergence angle is exactly as it is defined in the equation I posted, which is the same as in the figure Pietro posted. So if I understood correctly, it should be called "the half-divergence angle"? I can put that graph from the BBNBI paper to the documentation to make it clear.

[rui-coelho]

and why the factor 1/2 in the gaussian ? In Eq.1 from the paper of Otto it is absent....

[rui-coelho]

first test with new version shows that the bending is no longer present BUT.....to get the same deposition pattern as the ASCOT4/BBNBI4 run i need to divide the half-divergence by a further sqrt(2) see pictures below

BBNBI4

BBNBI5

[miekkasarki]

I had a typo in the equation I posted: it shouldn't be

$$e^{-\frac{1}{2}(\omega_x/\omega_v)^2}+e^{-\frac{1}{2}(\omega_y/\omega_h)^2},$$

but

$$e^{-\frac{1}{2}(\omega_x/\omega_v)^2-\frac{1}{2}(\omega_y/\omega_h)^2}.$$

As for the factor of $1/2$, my reasoning was that the above expression would then reduce to Otto's when $\omega_h = \omega_h$. But to be honest, I don't completely understood Otto's formula as integrating $\omega$ from $-\infty$ to $\infty$ does not yield 1...

However, if my reasoning was incorrect and the factor of $1/\sqrt{2}$ should be "inside" $\omega_{v,h}$ (so) then you Rui should be able to provide the half-divergence as is.

[miekkasarki]

So we just have to define what the exponential should look like and that defines our $\omega$.

[rui-coelho]

the point is really making sure the definition of half-divergence angle is consistent with the values we are given by the engineers. We need to ensure that they provide is:

1. Related to the variation over distance along beam path of the beam "radius"
2. This "radius" understood as the 1/e^2 folding of the intensity......or the power. This is a crucial point since if it is power then you square the expression and the exponential gets a factor 2 in the exponent....cancelling the 1/2 you have there...

Your probability gives intensity or power equivalent ?

[miekkasarki]

Now we have tests for the BBNBI and now I can say with confidence what the divergence is right now.

The top two plots show the vertical and horizontal divergencies. The blue curve is the fraction of particles at each angle, i.e. power (not intensity), and the black curve (and red) is the bi-Gaussian (halo+core) where both have the form:

$$\sim e^{-(\omega/\omega')^2}$$

$\omega'$ is the divergence that you use when providing the injector input so the value is not modified in any way within the code.

This means that the divergence that BBNBI5 uses:

• Corresponds to $1/e$ folding in power.
• It is therefore the "half-divergence angle" in engineering language.

Do you guys agree with me? If you want to verify it yourself, just run the BBNBI in this branch (I rebased this on the main so all the other fixes are now here as well).

[rui-coelho]

so, in theory i should only use the half-divergence angle and NOT apply an additional sqrt(2) ? I just raise this since it is a tricky point when we compare our NBI beams to usual laser beams. In laser beam optics the electric field goes like

E(x,y,z)~e^( - (x^2+y^2)/w(z)^2 ) w=beam "radius", easier to define once we define beam Intensity

Beam "intensity" is defined in W/m^2 and relates to the power in the electric field so

I(x,y,z)~e^( - 2*(x^2+y^2)/w(z)^2 )

So, the beam width (w) at any point of the propagation path (z) is defined as the radius at width the intensity reduces to 1/e^2 of the on-axis value. So, i am tempted to believe we need the 1/sqrt(2) in the divergence if we want to explain the matching curves you provided but please someone prove me wrong otherwise.....

@pietrovincenzi what do you think ?

[miekkasarki]

Yes right now as the code currently is, you shouldn't apply sqrt(2). Just give the half-divergence as is. Note that this was not the case earlier as I adjusted the definition in the code.