Simulations of Neutrino Factory Muon Capture with Chicane

Stephen Brooks (sb@stephenbrooks.org)   December 2003
1.  Introduction

The aim of a neutrino factory is to produce intense beams of neutrinos by trapping muons in a storage ring until they decay.  Muons, having a half-life τ½~1.56μs, are not readily available in nature and must be produced earlier on in the accelerator system.  The suggested way of doing this starts with a proton beam of several GeV and fires it at a target material, which emits pions (π+) that decay into the required muons (μ+) soon afterwards (τ½~26.0ns).

1.1.  Challenges

The main difficulty with this scheme is that the distribution of pions emitted from the source has an extremely large emittance: in the dataset used, 22% of the pions are actually going backwards in the lab frame (vz<0), so the emittance can hardly be calculated.  This means that any simulation based on conventional assumptions about the beam's general direction and size are likely to be inappropriate.  There are inevitably large losses associated with capturing this beam and focussing it into a form that could e.g. be accepted into a linac.

The decay of the pions into muons adds another unusual element to the system, as the other decay products are emitted at high speed and make significant changes to the original particle's mass and velocity.  The pion half-life corresponds to 7.8 light-metres, so the decays happen on the same timescale as it takes the particles to travel through the first few magnets.  This effect also makes simulation with a traditional tracking code potentially inaccurate.

1.2.  The RAL Design

This report focuses on the proposed design in which the pion source is a tantalum rod (20cm long, 1cm in radius) being hit by a 2.2GeV proton pulse, assumed to have a Gaussian longitudinal bunch shape of RMS one nanosecond.  The first component in the muon capture apparatus is a channel of transverse-focussing solenoids, the first of which encloses the rod (although the axis of the rod need not be parallel to the axis of the solenoid).  The lack of longitudinal focussing means that the particles with larger vz will work their way to the front of the bunch, which will have spread out somewhat by the end of the channel.  With this in mind, the next component is a 'chicane' of specially-designed bending magnets, in which the faster particles will go on paths further to the outside of two opposing bends in the accelerator.

Figure 1.  Schematic of the RAL muon front end design.
Figure 1.  Schematic of the RAL muon front end design.

The design is such that muons emitted directly down the solenoid channel with energies in range 120-270 MeV are meant to take an equal amount of time to get to the end of the bends.  This means that they are re-bunched longitudinally at that point, ready to be accepted into the third component of the design, an 88MHz linear accelerator (not yet simulated).

These components constitute a muon front end for the Neutrino Factory, with bunched muons leaving the system with a mean kinetic energy of 400MeV.  Other front end designs use different combinations of techniques to achieve the same output, for instance by using muon 'cooling' or several stages of linac acceleration.

2.  Computational Method
2.1.  Data Sources
  • Pion distribution emitted from the tantalum rod source, calculated from a Monté Carlo particle physics simulation of protons interacting with the rod [Paul Drumm].  This contains roughly 21000 CSV records, each specifying one pion's initial velocity and position at the rod surface.  Due to the nature of the simulation, these records also have 'weights' proportional to their relative probabilities, which are preserved throughout the rest of the simulation.
  • B-field produced by the bending magnets.  These magnets were purpose-designed to produce the desired variation in path-length as a function of muon energy [Mike Harold].  The field at any point within the magnets is not easy to calculate analytically, so instead the simulation interpolates a table of field values output directly from the magnet design package.  Due to the proximity of some magnets to each other, this data maps the fields through three adjacent magnets whose fringe fields interfere.
  • Placement of components along the beamline.  This is derived from a 'lattice file' specifying the structure of the accelerator, in combination with the list of optimisation parameters currently being tested (see §2.3).  Early versions of the program were less flexible with a hard-coded routine generating the beamline from the parameter list.  The new format shares the basic structure of inputs files to other codes (e.g. MAD), including features such as defining repeated blocks of components.  It also includes some less common features, such as notation to specify ranges of values to be optimised and (experimentally) a scripting language to deal with beamlines that depend on the parameters in a more complex way.
2.2.  Simulation

Space-charge effects are small in this system, due to the size of the beam (~10cm across) and low number of particles: five times fewer μ+ than there were p+ in the proton driver.  In the decay channel μ+ and μ are present in comparable numbers, further cancelling the effect, so the simulations did not include a space-charge solver.

2.2.1.  Basic Particle Tracking

This is done using a fourth-order Runge-Kutta (RK4) solver for each particle in the (x,y,z,vx,vy,vz) phase-space, all coordinates being relative to the lab frame.  A time-step was chosen by running batches of test particles through the solenoid channel and comparing the final positions with those obtained by halving the timestep.  With the 10-picosecond steps used in the real simulations, it is estimated that after 120 nanoseconds (approximately the total time the beam is simulated for), no errors in position of a particle will be more than 2.5mm, with RMS and mean errors 0.5mm and 0.3mm respectively.  This compares favourably with the size of the beam, which is typically at least 100mm across.

The way that forces are applied to the particles is special-relativistic and also fully integrated with the RK4 solver.  RK methods in general solve first-order differential equations, whereas equations of motion are second-order in normal space.  However, we find that the first time-derivative of the position in phase-space (d/dt)(x,v) = (v,dv/dt) is a function only of the original phase-space vector (x,v) - and possibly t, in the case of varying fields such as RF cavities.  The calculation of dv/dt is complicated as usual by the relativistic mass-increase, but eventually one can derive the formula dv/dt = (F−c-2(F·v)v)/m0γ, where F is the force as measured in the lab frame, i.e. q(E(x)+v×B(x)) with q=+e for our μ+ and π+.  These formulae are used in the program whenever the RK4 solver requests an evaluation of the derivative at some place in phase-space.

2.2.2.  Solenoid Fields

Currently the B-field contribution from solenoids is calculated using an analytic solution to Maxwell's equations with cylindrical symmetry, which yields the field as the sum of an infinite power-series.  An error analysis similar to that used for determining the RK4 timestep was done, which revealed that truncating this series at the cubic term gives maximum and RMS positional errors of 0.25mm and 0.01mm respectively, so including any higher terms would slow down the simulation with no significant gain in accuracy.

2.2.3.  Pion Decay Effects

The decay of pions is a random process in time, and the muon is emitted in an isotropically random direction in the rest frame of the original pion.  Both of these factors would be better suited to a regime in which probability-density functions fπ(x,v), fμ(x,v) of phase-space are simulated, rather than being sampled discretely by macro-particles.  While the development of such Vlasov solvers has been investigated, meaningful simulation of a 6D phase-space with reasonably high grid resolution has not been achieved.

To mitigate the noise introduced by the random decay of discrete particles, each macro-pion is decayed 10 times, producing 10 macro-muons of one-tenth the weight.  Simulations without this feature are moderately accurate but the optimiser often compares the merit of quite similar results, making random variations more of a problem.  When the bending chicane was added, the final muon transfer fell to the point where the optimisation was hindered by noise unless multiple decays were used.

2.2.4.  Delayed Emission

When calculating the longitudinal beam profile, there needs to be some mechanism to simulate that the original proton pulse is not instantaneous.  This is done by simulating several (usually 5) copies of the initial pions, but with the release of any one pion delayed by a random amount of time taken from the normal distribution t~N(μ=5ns,σ=1ns).  In cases where only the number of particles retained is important, all the particles may leave the source simultaneously, as without space-charge they do not interact.  (They may also be simulated individually).

2.2.5.  Absorbing Barriers

The simulation must take account of the various components situated in and around the beamline through which particles cannot pass.  One of these is the rod source itself; others are the regions outside solenoids and the rectangular apertures of the bending magnets.  Care must be taken that while stopping particles that pass outside a solenoid (for example) by a few inches, particles in a different part of the beamline intersecting the same plane (and thus technically 'outside' the solenoid) will not be removed.  If this is implemented wrongly, particles can make false collisions with components.

2.3.  Optimisation
2.3.1.  Recording Results

After each run of the simulation, the result can be appended to a text file in the two-line format shown below.

zsrc=584;rotsrc=080;d1l=322;s2f=997;s3f=883;s4f=577;snf=490;alt=415;narrow=658;na1=422;na2=962;na3=474;
2.060580 (20936 particles) [v4.11] {E416DD6E}

The first line contains a list of named parameters that determine the accelerator design within the range being considered in the current optimisation.  These have normalised values from 000 to 999, which can be rescaled onto the feasible range for that design parameter to find the actual value used.  They may also be used to switch between different schemes of designs - for instance, the 'narrow' parameter is used to switch between four different algorithms for varying the solenoid radii along the decay channel.

The second line starts with the 'score' for the accelerator.  From the optimiser's point of view, this just has to be a valid floating-point value with larger values corresponding to better designs (e.g. percentage efficiency).  The rest of the line shows the number of particles used, the version number of the simulation software and a checksum of the beamline design file (a comment may also be added).  This information is useful for distinguishing relevant results from those generated under different conditions or earlier (perhaps incorrect) versions of the software.

2.3.2.  The Algorithm

Optimisers are often implemented with steepest descent or hill-climbing algorithms that use the derivative of the score with respect to the parameters to determine which direction the optimiser will search along next.  Unfortunately effects from pion decays add a small random variation to the score at any fixed design, making the calculation of derivatives impractical.

Instead, the optimiser uses an idea based on genetic algorithms (GAs).  A conventional GA starts by running a number of random designs (the first generation) and produces a new generation using combinations of features from the best in the previous one, plus a small amount of random 'mutation'.  This process is repeated until the score no longer increases.  That would be far too wasteful for the purposes of this simulation, however, as it throws away the results from each generation after it generates the next (and a single simulation takes several minutes).  Instead the optimiser logs all of the designs simulated (good or bad) to a file, and chooses a new design based on that data.

When the number of records in the file is small (e.g. less than 10), the program simulates a random design.  Otherwise, it chooses one of the following four options with equal probability.

  • Mutation.  Simulate a design that is similar to a previous good design, but has some parameters varied by random amounts.  This is useful for fine-tuning a near-optimal design.
  • Crossover.  Pick two good previous designs and interlace the parameters from both, so the resulting design has each value chosen randomly from one of the two originals.  Crossover has the potential to 'combine the best aspects' from two different designs.
  • Linear interpolation.  Generate a design that lies somewhere on the straight line between two previous good designs in the 'optimisation space'.  The other two methods vary the parameters quite separately, but this will vary them all in a correlated way.
  • Random.  Random designs are still simulated 25% of the time to ensure that the optimiser does not 'miss' a global minimum by getting stuck in a confined region with a local minimum.
2.3.3.  Distributed Computing

Since the results from independent optimisations are all still valid, this optimiser can be deployed as a distributed computing project, with near-independent calculations being made on many different machines.  This resulting project status page is currently at http://stephenbrooks.org/muon1.  At times a "best results" or "results sample" file for the whole project is made available so that those who download it can base their optimisation on results found by other computers.

3.  Optimisation of the Decay Channel
3.1.  Original Design

Before work on these simulations started, the muon front end had been developed [Grahame Rees] using standard modelling codes such as TRACE3D and Parmila.  The basic design was as follows.

Element
S = Solenoid, D = Drift space
Length (m)Field (T)Aperture Radius (m)
S10.406620.00.1
D10.5718
S20.4-3.30.3
D20.5
S30.44.00.3
D30.5
S40.4-3.30.3
D40.5  
S50.43.30.3
D50.5  
S6-D239 more repeats of the above block
S240.4-3.30.3
D240.5
S250.43.30.15
D250.5  
S260.4-3.30.15
D260.5  
S27-D344 more repeats of the above block

The tantalum rod pion source is positioned concentric with S1 and rotated by 0.1 radians about the Y-axis.  This rotation reduces the number of pions lost from reimpacting the rod after one turn of a helical path in the solenoid field.  Some pions are moving forward rapidly enough to clear the rod before returning to the same transverse location, but for those that are not, the shifting transverse footprint of a tilted rod with distance makes it possible to avoid a collision.  Having the rod at too small an angle is also undesirable when wanting to remove the remaining proton beam, which has not interacted with the target, from the system.

The results from one preliminary simulation to determine the effect of rod angle on pion capture are graphed in figure 2.

Figure 2.  How rod angle affects the performance of S1.
Figure 2.  How rod angle affects the performance of S1.

This simulation was made with only the solenoid S1 in place, so the 'beam kept' figure counts the particles emerging from the correct end of this solenoid and 'wrong way' counts those from the opposite end.  The graph suggests that a rod angle of ~0.45 radians would produce the lowest loss and avoid high re-interaction between the emitted particles and the rod.

For reference, the original design was tracked with the code 36 times and yielded a transmission efficiency of 3.11±0.01% μ+ out per π+ emitted from the rod, where the range quoted is a 95% confidence interval over the variation caused by random simulated pion decays.  Conversions to other figures of merit for all efficiencies quoted here are found in the appendix, while units of μ+/π+ will be used here.

3.2.  Parameters Varied

Two methods of optimising the solenoid channel have been investigated.  The first generalises the design given in §3.1 by varying the 12 parameters shown below.

Name(s)DescriptionRange of valuesNotes
zsrcZ-displacement of the rod source-0.15 to +0.15 metresMeasured as the position of the centre of the rod relative to the centre of S1.
rotsrcRotation of the rod source about the Y axis0 to 0.5 radians
d1lLength of drift space D10.5 to 0.6 metresThe drift spaces in this area of the accelerator are primarily determined by other instrumentation, but this one can be varied slightly.
s2f s3f s4fThe strength of solenoids S2, S3 and S4±2.5 to ±5.0 TeslaThe sign as well as the strength can be changed by the optimisation; very low fields are not worth considering.
snfMagnitude of field in solenoids S5-S34±2.5 to ±4.0 TeslaS5's field has the same sign as this.  The fields in later solenoids will have the same magnitude but sign controlled by alt.
altControls how the solenoid fields S6-S34 change signDiscrete; two valuesEither the solenoids alternate in sign (+-+-+-), or they alternate in groups of five (+++++-----+++++-----).
narrowMethod of narrowing the solenoids S10-S30Discrete; four valuesThe solenoids narrow from 30cm to 15cm radius in one step; two steps via 22.5cm; three steps via 25cm, 20cm; or continuously in a section of linearly-decreasing radii.
na1 na2 na3Solenoid narrowing pointsSolenoid number; integer from 10 to 30These control the points at which the step-narrowings take place, or where the linear narrowing begins and ends.  na2 and na3 may be unused.

When this phase proved successful, the design range was extended to ~130 parameters by varying every component length, field and radius independently.  The following practical restrictions defined the design space:

  • Solenoids can be from 20 to 60cm long, apart from S1 whose upper limit is 45cm.
  • Drift lengths range from 50cm to 1 metre.
  • The total length of the lattice is kept close to 30 metres by asserting that the final solenoid is the first whose start is more than 30 metres from the start of S1.
  • The radius of the superconducting S1 is 10cm and the radius of the final solenoid must be 15cm.  The rest can have radii from 10 to 40cm.
  • Fields strengths of either sign up to 4 Tesla are allowed for most solenoids and 5 Tesla for S2-S4.  S1 has a positive field of between 0 and 20 Tesla.
  • The rod can be placed in roughly the same range of locations as before.

The first optimisation is a subset of this range, so the second ought to produce a yield at least as good, with the amount of improvement indicating how suitable the extra design assumptions in the first phase were.  The limiting case is when no improvement is made, meaning the optimal design lies within the first subset.

3.3.  Results from 12-Parameter Optimisation

The highest-scoring result recorded was the following.

Parameter(s)Value(s)
zsrc rotsrc-0.01006 metres, 0.03253 radians
d1l0.5172 metres
s2f s3f s4f+4.970, +4.660, +3.849 Tesla
snf alt+3.985 Tesla, signs alternate in groups of five solenoids (+++++-----+++++ etc.)
narrow na1 na2 na3Solenoid radii decrease in three steps: to 25cm at S10, 20cm at S13 and to 15cm at S30.

This was given a score of 7.066% μ+/π+, but due to statistical effects (and the fact this was the best result, so likely having an upward fluctuation), the true figure is probably closer to 6.5%.  These effects were more severe than in the 'high-generality' optimisations discussed later because an earlier version of the code, without multiple decays for each pion, was being used.

Unlike the preliminary simulations, the optimisation of the full decay channel seemed to favour small rod rotation angles, perhaps because capture of the pions by S2 and S3 is more difficult when the rotated rod produces a larger X-axis spread of particles.  The tendency towards choosing the strongest possible focussing solenoids for the majority of the beamline is not surprising given that the goal is to minimise beam loss, however the low value of s4f suggests a higher field at that point might not be beneficial.

3.4.  High-Generality Optimisation

Moving onto a larger design space, the program reached a design whose transmission was confirmed by multiple runs to be 9.65±0.01% μ+/π+.  As before, many of the parameters were near the ends of their allowed ranges.  However, not all were exactly at the extreme value, which is understandable given how difficult a yield change caused by a <0.01T field variation in one magnet is for the optimiser to detect in the presence of noise.  Thus to simplify the design presented here, parameters within 10 quanta of the ends of the range were rounded to the ends.  The resulting design also scored 9.65% when tracked and is listed below.

ElementLength (m)Field (T)Aperture Radius (m)
S10.4520.0[0.1]
D10.5
S20.65.00.3514
D20.5
S30.42305.00.4
D30.5
S40.38064.1890.4
D40.6612
S50.52994.00.4
D50.5075
S60.63.8240.4
D60.5170
S70.64.00.4
D70.5  
S8-D2518 more repeats of the above block, apart from three exceptions:
S130.59564.00.4
S140.64.00.3862
D200.5175  
S260.63.9040.4
D260.6041
S270.64.00.4
D270.5
S280.64.00.3159
D280.5
S290.24.0[0.15]
Parameters not at the extremes of their range are shown in italics.  Those fixed by the optimisation are shown in [brackets].

The tantalum rod in this design is rotated by 0.09359 radians, with centre 0.2252 metres from the beginning of S1, which makes it nearly concentric as S1 is 0.45 metres long.

Figure 3.  Shape of the pion beam on entering S2.
Figure 3.  Shape of the pion beam on entering S2.

As before, the optimiser has chosen to have the strongest possible focussing for the majority of the beamline by using the longest solenoids with the maximum field and minimising the length of the drifts.  The first few components have been optimised in other ways, perhaps to account for the unusual shape of the beam in the early stages of the channel (see figure 3).

The strength of the field in S4 is substantially lower than the maximum allowed, just as happened with the choice of s4f in the last optimisation, which indicates this is a genuine feature and not just the optimiser converging badly.  As another check, the region from S5 to D28 in the above design was 'homogenised' to look entirely like the repeating block and tested with the code.  This produced a yield of only 9.09%, which is significantly below the original.  However, when the three 'exceptions' from the middle of the beamline were removed, a short re-test produced 9.68%, meaning that they could probably be omitted.

Perhaps the most significant change of the lattice has been in the polarity of the solenoids: while the previous optimisation had the sign change every five solenoids instead of the alternation present in the original, this one has taken the process to its conclusion and does not have any changes of sign.  Considering the transverse motion of particles, we see that particles travel in straight lines in drift spaces and around circular arcs (in the transverse projection) while in solenoids.  For a fixed energy, there are closed transverse paths possible for both the alternating (FODO) and the same-sign (FOFO) solenoidal lattices, resembling a bow-tie and an oval racetrack respectively.  Higher energy particles will have a larger radius of curvature in the solenoids and will turn through a smaller total angle (and lower energy particles the opposite).  This means that the path of an off-energy particle will not close.  In the FOFO lattice this causes the racetrack orbit to precess because of the extra rotation on each circuit.  However, due to the opposing bends in the FODO bow-tie, this orbit suffers a translation each time it fails to close, shifting off-energy particles ever more off-axis and dispersing the beam.  The precession process, by contrast, remains bounded in transverse space since any translation that may be present in the orbit defect will eventually be cancelled when it has rotated by 180°.
The way this was found illustrates the ability of a high-dimensional optimiser to generate solutions qualitatively different to the reference design, sometimes even incorporating design principles that (as in this case) were not foreseen by the authors.

4.  Simulation with the Bending Chicane

The second optimisation presented here varies the solenoid channel in the same way as before but also tracks the particles through the second segment of the design.  The score for each configuration is calculated from the number of particles retained at the end of this second part, meaning that although the chicane itself cannot be modified due to beam optics considerations, the decay channel can be refined to be optimally compatible with it.

4.1.  Structure of the Bend

The bend consists of four identical 'triplets' of magnets (B1, B2, B1) with fields such that a pencil beam of muons with varying energies emerges from a triplet as a dispersed beam and rotated by 101°.  When that dispersed beam enters a second triplet at the correct angle, it is turned back into a pencil beam (up to the correctness of the field within the magnets) but by then the higher-energy muons will have travelled further, giving the slower ones time to catch up.  To cancel the dispersion produced by the decay channel, 404° of bends are needed, so to avoid layout problems the second pair of triplets bends in the opposite direction.

Figure 4, produced using the same tracking code as used for the optimisation, shows how the magnets will bend the paths of muons with energies of 120 MeV (red) up to 270 MeV (violet) in steps of 10 MeV.

Figure 4.  The reverse phase slip bending chicane in action.
Figure 4.  The reverse phase slip bending chicane in action.

Not all the particles make it through the chicane: the iron poles of the magnets are at roughly ±15cm from the mid-plane, close enough to create losses, especially since vertical focussing is quite weak.  This is modelled by placing aperture restrictions at either end of each sector magnet.

The reference design for the solenoid channel (§3.1) was tested with the chicane added and scored 1.127±0.006% μ+/π+, giving the chicane an apparent efficiency of 36%.

4.2.  Longitudinal Distribution

A preliminary simulation calculated the distribution of muon arrival times at the end of the bending chicane, to see if the beam will be sufficiently bunched to be accepted into the 88MHz linac.

Figure 5.  Longitudinal beam profile upon re-bunching by the chicane.
Figure 5.  Longitudinal beam profile upon re-bunching by the chicane.

As figure 5 shows, the bunch emerging from the chicane is approximately the same length as the linac bucket (an 88MHz wavelength is shown for comparison), though some loss of the distribution tails will occur.
A lot of the apparent bunch length reduction in the chicane is due to energy selection by the horizontal aperture of the magnets when the beam is dispersed, since the particles with extreme energies tend to be at the ends of the bunch.  For example, in one configuration the raw output of the solenoid channel has an RMS duration of 9.17ns, which reduces to 2.88ns for the 59.9% of the beam that is within the chicane's design acceptance of 120 to 260 MeV.  If only the sub-120 MeV stragglers are removed, the figure is just 3.40ns for 90.7% of the beam.
On leaving the chicane, the bunch has an RMS duration of 2.11ns in which 93.1% of the particles are of the design energy.  If only particles in that range are considered, the bunch extent has reduced by 30% from 2.88 to 2.01ns, which is due mainly to the dynamics of the chicane rather than losses.

4.3.  Results from 12-Parameter Optimisation

The design from the 12-parameter optimisation of the decay channel (§3.3) was run with the bending chicane appended, to provide a comparison between separate and joint optimisation, and produced a yield of 1.878±0.027% μ+/π+.  Then the full optimisation, scoring by transmission through the end of the chicane instead of just the solenoid channel, generated a design with a transmission of 1.930±0.013%.  This is higher by a statistically significant amount, but the code without multiple pion decays was still being used at this point, and with the lower yields coming from the chicane, noise started to impair the optimiser, so these figures may not be quite optimal.

4.4.  High-Generality Optimisation

Increasing the number of parameters improved the transmission to 2.704±0.004% μ+/π+.  One comparable result for separate optimisation is the design of §3.4 with the chicane added, which scores only 2.004±0.005%.  Closer examination reveals that the restriction of fixing the final solenoid's radius in this lattice (which was meant to be a substitute for the finite aperture of the chicane) actually produces a betatron focus there when optimised, in which the transverse velocities of the particles are still large.  Removing this restriction gives a decay channel transmitting 15.80±0.01% on its own and 2.238±0.004% through the chicane.  Still, jointly optimising the two sections has produced an improvement of just over 20%: not insignificant given that the full neutrino factory has many sections and a compounded 20% improvement at each interface would correspond to an overall factor of 2× or 3×.  This performance would only be accessible via holistic approaches such as optimising using a (currently hypothetical) code capable of tracking through the whole machine, or tuning the magnet strengths on the real machine after it is built, the latter being able to recover only the fraction not requiring machine geometry changes.

Only a few of the components in the new design deviate from the pattern of using the highest-possible solenoid fields and the shortest-possible drifts for maximum average field on-axis.  There are some weaker solenoids near both ends of the channel, as has happened before for matching with the input and output.  The only non-minimal drift could well be an artefact of the optimiser having not quite converged.

SolenoidField (T)DriftLength (m)
S120D10.5
S25D20.5
S33.088D30.5
S44.429D40.5
S5-S254D5-D250.5
S263.752D260.5
S274D270.5345
S283.760D280.5
S29-S314D29-D310.5
S323.431D320.5

The other parameters in this design are of more interest.  In the decay-channel-only simulations, the program predominantly opted for the (logical) approach of using the largest possible solenoid radii.  By contrast, the new design has an evident solenoid radius structure (see figure 6).

Figure 6.  Solenoid radii in the high-generality optimisations.
Figure 6.  Solenoid radii in the high-generality optimisations.

The decreasing trend is fit well by the linear formula 0.402−0.003n for the radius of Sn from S7 to S28.  Note that the optimisation algorithm itself does not treat the parameters as being in any particular order, so trends such as this cannot be artificial and must come from the beam dynamics.  Secondly, note that if hard-ended solenoids had been used in the simulation, narrowing the radii in this way would have produced a result no better than keeping them at maximum size until the very last solenoid, because particles would experience the same fields in such a linear model in either case until they are lost.  This suggests that the nonlinear part of the solenoid end-fields is being used to achieve some sort of additional focussing.
The fact that the solenoids have not narrowed all the way to 15cm is another indication that the use of that design constraint in §3.4 was not appropriate.

Figure 7.  Solenoid lengths in the high-generality optimisations.
Figure 7.  Solenoid lengths in the high-generality optimisations.

Figure 7 shows that the optimiser has settled on an intermediate value for the solenoid lengths in the main part of the channel, in contrast to how it made them all as long as possible without the chicane.  The median of the lengths used is 0.463m.
Some analysis of regular solenoid focussing channels using linear optics was done in order to explain this figure.  The Larmor rotation within each solenoid is given by

θ =
qBzL
m0γvz
,q, m0 = Particle charge, rest mass
Bz, L = Solenoid field, length
γ, vz depend on particle velocity

which equals 160° in the case of 120MeV muons down to 90° in the case of 260MeV muons.  Angles that are too small (or more generally, too close to a multiple of 360°) will mean particles will only have rotated slightly each time the come to a drift space, leading to a large transverse orbit.  The FOFO lattice is also not unconditionally stable, as if the drift length ratio D:= Ldrift / Lsolenoid is greater than a certain value

Dcritical =
4
θ
 
sin θ
1−cos θ
+
2
1−cos θ
 
,

the transfer matrix becomes unbounded when iterated.  For this lattice, D=1.08 so we must have θ<192° for stability (discounting solutions θ>360°).  The range 90°-160° can be seen as avoiding both of these problems by a suitable margin.

Finally, to check how much of the optimiser output is "noise" surrounding the trends already identified, rounded versions of the lattice were checked, described by the table below.

SolenoidLength (m)Field (T)Aperture Radius (m)
S10.394920[0.1]
S20.650.3093
S30.52353.0880.402-0.003n
S40.4634.429
S54
S610.463, 20.5584
S7-S300.463
S310.3461
S320.46513.431

All drifts were set to 0.5m long.  There was some doubt as to whether the particularly long length of S6 (see figure 6) was genuine or an optimiser error.  In the first simulation, the solenoid lengths were rounded as far back as S4 but this decreased the transmission to 2.678±0.008%, which is lower than the original with a statistical significance of 5.7σ.  In the second rounded design the value of S6 (but not S5 or S4) was restored, which gave an efficiency of 2.704±0.006%, just as good as the non-rounded design.

5.  Future Work

Recently interest has moved onto a new sequence of components after the decay channel [Grahame Rees], consisting of a muon cooling ring with absorbers to reduce the emittance of the beam over several turns, preceded by an RF system to phase-rotate the decay channel output into a smaller energy spread.  At the other end of the beamline, it has become clear that various practicalities demand the target area be changed from the model used here.  These and other work to be completed with this code are listed in detail below.

5.1.  Extension to RAL Muon Cooling Ring
  • Add RF cavities to the simulation: these are needed in both the cooling ring lattice and also in the remaining 88MHz muon linac of the chicane lattice.
  • Include a model of muon scattering in the cooling ring absorbers (which will probably be liquid H2).  Both constant-thickness and unusual wedge-shaped absorbers are required for 6D phase-space reduction.
  • The RF phases, absorber thicknesses and shapes should all be optimised along with the usual parameters for the solenoids and drifts in the rest of the lattice.
5.2.  Developing a Consistent Target Area
  • The first solenoid (S1) should be enlarged to allow more pions to rotate within it and to reduce heating by losses on its inner surface.
  • S1 may also have to be split into two solenoids in a 'Helmholtz coil' arrangement in order for targets on a rotating wheel to enter from the side (a single tantalum rod cannot dissipate more than 30kW of heat, we require at least 1MW).
  • The rod is not a stopping target and so several megawatts of unused protons must be taken to a beam dump.  Consideration must be given to the geometry involved and how this beam can be stopped from hitting one of the solenoids.
  • Repeat the simulation using the pion distribution from an 8GeV instead of a 2.2GeV p+ beam, to account for a more likely design of the RAL proton driver and to explore the sensitivity of the optimal lattice to the initial pion distribution.
5.3.  Modifications to the Simulation Code
  • Use variable-length timesteps if possible, to increase the speed or accuracy of the particle tracking.
  • Check to see if space-charge effects become significant when the muon bunches are cooled to a higher density.
  • Consider solenoids as having thickness rather than being perfect cylinders.
5.4.  Long-Term Goals
  • When the code is developed enough, run optimisations on several proposed front ends (from FNAL, CERN and others) to evaluate the performance of contrasting designs, some with and some without muon cooling.
  • Investigate other exotic components such as magnetic horns.
A.  Conversion of Transfer Efficiencies

The simulations described here used a dataset generated for 105 2.2GeV protons hitting the rod.  The total π+ weight produced was 21285.3, so conversion is done at 0.212853 π+/p+.

BeamlineEfficiency in μ+/π+
95% C.I. uncertainties for last digits in brackets
μ+ per proton on rodμ+ per proton.GeV
Original decay channel3.114(10)%0.6627(22)%3.012(10)×10-3
12-parameter optimised~6.5%
Simulated with earlier, less-accurate version of the code
~1.4%~6.3×10-3
High-generality optimised9.645(13)%2.053(3)%9.332(12)×10-3
Original to end of chicane1.127(6)%0.2399(13)%1.090(6)×10-3
12-parameter optimised1.930(13)%0.4107(28)%1.867(13)×10-3
12-parameter for decay channel through chicane1.878(27)%0.3998(57)%1.817(26)×10-3
High-generality optimised2.704(4)%0.5756(8)%2.616(4)×10-3
High-generality for decay channel through chicane2.004(5)%0.4267(11)%1.939(5)×10-3