11institutetext: Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA

# Protection Related to High-power Targets

M.A. Plum

###### Abstract

Target protection is an important part of machine protection. The beam power
in high-intensity accelerators is high enough that a single wayward pulse can
cause serious damage. Today’s high-power targets operate at the limit of
available technology, and are designed for a very narrow range of beam
parameters. If the beam pulse is too far off centre, or if the beam size is
not correct, or if the beam density is too high, the target can be seriously
damaged. We will start with a brief introduction to high-power targets and
then move to a discussion of what can go wrong, and what are the risks. Next
we will discuss how to control the beam-related risk, followed by examples
from a few different accelerator facilities. We will finish with a detailed
example of the Oak Ridge Spallation Neutron Source target tune up and target
protection.

Keywords
Target protection; beam; proton; high-intensity.

##  0.1 Introduction to high-power targets

In general, high-power targets are needed to create either lots of secondary
particles for applications where the interaction cross-sections are low, or to
make as many secondary particles as possible when the creation cross-sections
are low. Example applications include:

  * •

neutron spallation targets (LANSCE, ISIS, SNS, J-PARC, PSI);

  * •

muon production targets (J-PARC, PSI, TRIUMF, ISIS);

  * •

Isotope Separation On-Line (ISOL) facilities (CERN, TRIUMF, IRIS);

  * •

material irradiation studies (IFMIF);

  * •

antiproton production (FNAL);

  * •

neutrino production (FNAL, J-PARC, CERN).

High-power targets come in many shapes and sizes, and have many uses. Some
example high-power targets are shown in Figs. 1–3, and parameters [1] from
some example high-power facilities are shown in Table 1. High-power targetry
is a highly developed and complex field, and there are many technological
challenges. Some of the biggest beam-related challenges are:

  * •

removing the heat generated mainly by the beam, but also by nuclear decay;

  * •

mechanical shock due to thermal stress from pulsed beams;

  * •

radiation damage, including swelling and embrittlement;

  * •

target handling (including installation and maintenance in a high-radiation
environment);

  * •

corrosive environment;

  * •

beam parameters—matching target requirements to what the accelerator can
deliver.

Figure 1:  Some example neutron spallation targets. Images reproduced from Refs. [6, 2, 3, 4, 5] Figure 2:  Some example neutrino and pbar targets. Images reproduced from Refs. [8, 7] Figure 3:  Some example isotope production targets. Images reproduced from Ref. [9] Table 1: Some example high-power target facilities Facility |  Status |  Target material |  Beam pulse duration |  Rep rate |  Proton energy |  Time avg beam power |  Peak time ave power density |  Peak energy density
---|---|---|---|---|---|---|---|---
|  |  |  ($\mu$s) |  (Hz) |  (GeV) |  (MW) |  (GW/m3) |  (MJ/m3)
ISIS |  Operating |  W |  0.4 |  50 |  0.8 |  0.16 |  0.25 |  5
LANSCE–Lujan |  Operating |  W |  0.3 |  20 |  0.8 |  0.16 |  0.5 |  25
NuMI |  Operating |  C |  8.6 |  0.53 |  120 |  0.4 |  0.32 |  600
SINQ/Solid target |  Operating |  Pb–SS clad |  CW |  |  0.57 |  1 |  1 |  NA
SINQ/MEGAPIE |  Completed |  Pb–Bi |  CW |  |  0.57 |  1 |  1 |  NA
JSNS |  Operating |  Hg |  |  25 |  3 |  1 |  0.63 |  25
SNS |  Operating |  Hg |  0.7 |  60 |  1 |  2 |  0.8 |  13
ESS—long pulse |  Proposed |  Hg |  2000 |  16.7 |  1.3 |  5 |  2.5 |  150
ESS—short pulse |  Proposed |  Hg |  1.2 |  50 |  1.3 |  5 |  2.5 |  50
EURISOL |  Proposed |  Hg |  3 |  50 |  2.2 |  4 |  100 |  2000
IFMIF |  Proposed |  Li |  CW |  |  0.04 (D2) |  10 |  100 |  NA
LANSCE–MTS |  Proposed |  Pb–Bi/W |  1000 |  120 |  0.8 |  0.8 |  2.4 |  20
US Neutrino Factory |  Proposed |  Hg or C |  0.003 |  15 |  24 |  1 |  3.8 |  1080
AUSTRON |  Proposed |  W |  1 |  10 |  1.6 |  0.5 |  |

##  0.2 Target environment

The target environment is often rough vacuum or helium. Vacuum is good because
it does not interfere with the beams (primary or secondary). Helium atmosphere
is good because it helps cool components while minimizing beam scattering;
however, impurities in helium can lead to corrosion of components. Beam
transport lines leading up to the target need good vacuum to minimize beam
loss due to scattering, so high-power targets often involve windows to
separate the beam line vacuum from the target environment. Safety separation
must also be considered (usually required for a liquid metal target). A
partial pressure of air in the target environment can result in the formation
of nitric acid, which can cause stress corrosion cracking in high-strength
steel. Nitric acid can also cause vacuum leaks in thin bellows exposed to
stray beam with air outside. For example, IPNS had a target bolt fail from
stress corrosion cracking in a nominal helium atmosphere with air impurities
when high-strength steel was substituted for stainless steel. Also, ISIS
observed corrosion around the target/reflector assemblies and consequently
limited impurities to below a couple of percent in helium. And, as shown in
Fig. 4, Fermilab had a broken chain due to acidic vapour/condensate from air
ionization.

Figure 4:  High-strength broken steel chain from hydrogen embrittlement caused
by acidic vapour/condensate from air ionization (MiniBooNE 25 m absorber).
Figure reproduced from Ref. [10].

##  0.3 Beam windows

Just like the target, beam windows are also challenged by high-power beams,
due to heating, thermal stress, radiation damage, etc. An example beam window
is shown in Fig. 5. This window is used to separate the SNS beam line vacuum
of approximately $1\times 10^{-7}$ Torr from the target environment, which is
1 atm He. The window is made of water-cooled Inconel.

|
---|---
Figure 5:  Side view and end view of the SNS proton beam window, located about
2 m upstream of the neutron spallation target. The end view also shows the
eight thermocouple halo monitors; the red ellipse indicates the approximate
beam size. Images reproduced from Refs. [4, 11].

##  0.4 Beam parameters

The highest-power targets push the edge of achievable technology, and to live
on the edge requires a very narrow range of beam parameters to avoid
overpowering the target. High-power targets require the nominal beam
distribution to be in the correct location, without exceeding the maximum beam
current or the maximum beam density. The most important beam parameters are:

  * •

beam position;

  * •

beam size and shape (i.e. distribution);

  * •

beam energy;

  * •

beam current;

  * •

beam pulse length, repetition rate, energy per pulse (short pulses cause a
pressure pulse).

High-power targets often prefer flat or uniform beam distributions because
they minimize the beam density and thereby minimize the density of energy
deposited in the target. However, unless special measures are taken in the
accelerator/beam delivery systems, the beam distribution will usually be
nearly Gaussian. A possible target protection requirement could then be to
monitor the beam distribution and to shut off the beam if the distribution
exceeds requirements. For example, Fig. 6 shows two possible beam
distributions. One is rectangular and the other is Gaussian. They both have
the same area (or beam power), but the Gaussian distribution has twice the
peak beam density.

Figure 6:  Two possible beam distributions. One is rectangular and the other
is Gaussian. They both have the same area (or beam power), but the Gaussian
distribution has twice the beam density.

###  0.4.1 Methods to flatten the beam distribution

Depending on the application, there are several ways to flatten the beam
distribution, to make it more favourable for the target by decreasing the peak
density. For storage rings and synchrotrons that employ multiturn injection
(e.g. LANSCE, SNS, J-PARC), injection painting can be an easy way to control
the distribution. Another method is to employ multipole magnets (e.g.
octupoles) in the beam transport line to target. This method was recently
implemented at J-PARC, as shown in Fig. 7. Note that the beam must be well
centred in the octupoles for proper flattening, since off-centre beams will
produce skewed distributions. Another method is a rastering system, that
‘draws’ the desired profile by quickly moving a small beam spot on the surface
of the target. This is the method planned for the ESS, as shown in Fig. 8. All
these methods rely on magnets or pulsed power supplies to achieve the desired
beam profile. This highlights possible equipment that would be important to
include as inputs to the machine protection system (MPS), and also possible
beam monitoring systems, such as beam distribution monitors, that could also
be part of the MPS interlock system.

|  |
---|---|---
Figure 7:  Simulations of beam distributions on the J-PARC neutron production
target. Left: no octupoles. Centre: with octupoles. Right: with octupoles and
beam off centre at octupoles. Figures reproduced from Refs. [12, 13].
---
Figure 8:  The ESS design to achieve the desired beam distribution on the
target by rastering the beam. Top: the beam transport line just upstream of
the target. Bottom: simulation of the achieved beam distribution. Figures
reproduced from Ref. [14].

##  0.5 What can go wrong?

In this section, we briefly describe some examples that highlight the
importance of MPSs for high-power target protection.

In May 2013, a gold target (shown in Fig. 9) melted in the Hadron Experimental
Facility at J-PARC and radioactive material was released from the building
[15]. The root cause was traced to a quadrupole magnet power supply in the
main ring, which caused the beam to be spilled in 5 ms instead of the nominal
2 s. Due to this accident all beam operations were terminated for about 8
months and beam operations to the hadron hall were stopped for more than one
year.

In January 2010, the neutron production target (shown in Fig. 10) failed at
the ISIS second target station [16]. The root cause was traced to a
combination of the high-density beam profile and a water leak. The water
disassociated in the high-temperature environment and then the oxygen attacked
the grain boundaries on the tantalum cladding. The cladding essentially fell
apart, leaving a hole in the target. In 2004 the beam density on the PSI
neutron production target increased to 3.5 times the nominal value [18, 17]
because the transport line (shown in Fig. 11) was accidentally set up for the
case of the muon target being inserted, but actually the muon target was not
inserted. The over-density condition was caught by a newly installed beam
distribution monitor (VIMOS) before any serious damage could occur.

We also need to be concerned with beam windows. They are often part of the
target system, and can also be damaged by the beam. High-power beam windows
must be cooled, usually with water. Windows can fail due to over-focused or
off-centre beams. Radioactive water can be spilled into the beam line vacuum
and/or the target environment. In 1996 the beam window at the LAMPF linac beam
dump failed [19], probably due to a combination of an unusually small beam
size (set up on purpose for an experiment) and a weld joint where a
thermocouple was attached to the air-side surface of the water-cooled window.

Figure 9:  Photographs of the gold target in the J-PARC hadron hall, melted by
a proton beam that was accidentally extracted over 5 ms. Figure reproduced
from Ref. [15]. Figure 10:  A photograph of a failed target from the ISIS
second target station. The tantalum cladding on the tungsten target failed.
Figure reproduced from Ref. [16]. Figure 11:  The beam transport line leading
to the PSI neutron production target, showing the beam sizes for the two
different sets of parameters (one for the muon target inserted, the other for
the muon target retracted). Figure reproduced from Ref. [17].

##  0.6 Calculating the beam density

To demonstrate how the beam density depends on various beam parameters, we
will consider several examples.

###  0.6.1 Example 1: rectangular DC beam

Assume a 2 MW DC beam of 1 GeV protons with rectangular cross-section 2 cm
$\times$ 2 cm. Assume that it stops 0.75 m into the target. Also assume
uniform energy deposition and no scattering. The energy deposition will be as
depicted in Fig. 12. The 1-s power density is then

$(2\,\mathrm{MW})/(2\,\mathrm{cm})/(2\,\mathrm{cm})/(0.75\,\mathrm{m})=6.67\,\mathrm{GW/m}^{3}.$

Now assume that if the beam is too far off centre it must be shut off before
it can deposit 100 J. What is the acceptable time delay (or MPS response
time)? The beam must be shut off within (100 J)/(2 MJ/s) = 50 $\mu$s.

Figure 12:  A beam pulse with a rectangular cross-section

###  0.6.2 Example 2: rectangular pulsed beam

Consider a pulsed version of the DC beam in Example 1. Assume that the beam is
divided into 10 pulses per second, and the duration of each pulse is 1.5
$\mu$s.

  * •

The 1-s power density is still (2 MW)/(2 cm)/(2 cm)/(0.75 m) = 6.67 GW/m3.

  * •

The energy per pulse = (2 MW)/(10 Hz) = 0.2 MJ.

  * •

The energy density per pulse in target = (0.2 MJ)/(2 cm)/(2 cm)/(0.75 m) = 667
MJ/m3.

  * •

The peak power density = (667 MJ/m3)/(1.5 $\mu$s) = $4.45\times 10^{14}$ W/m3.

Each beam pulse contains more than 100 J of energy, so even one off-centre
pulse will exceed the 100 J limit in Example 1. With 1.5 $\mu$s beam pulses,
it is not practical to shut off the beam mid-pulse. We can only prevent the
next pulse from occurring, and there is a fair amount of time to do that—about
0.1 s.

###  0.6.3 Example 3: pulsed Gaussian beam

Consider a Gaussian version of the pulsed beam in Example 2. Assume that the
Gaussian shape is characterized by $\sigma_{x}=\sigma_{y}=0.5$ cm, as shown in
Fig. 13. The functional form of the Gaussian is

$f(x,y)=A{\rm
exp}\left({-\left(\frac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}+\frac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}\right)}\right).$

The volume under this two-dimensional Gaussian is

$V=\int_{-\infty}^{\infty}{\rm d}x\int_{-\infty}^{\infty}{\rm d}yf(x,y)=2\pi
A\sigma_{x}\sigma_{y}$

and the 1-s peak power density is

$(2\,\mathrm{MW})/(2\pi)/(0.5\,\mathrm{cm})/(0.5\,\mathrm{cm})/(0.75\,\mathrm{m})=17\,\mathrm{GW/m}^{3},$

which is 2.25 times greater than the rectangular distribution case. The per-
pulse peak energy density and peak power density also increase by a factor of
2.5.

Figure 13:  A two-dimensional Gaussian beam profile

##  0.7 Target protection

We will consider only the beam-related protection. The targets themselves will
have their own separate non-beam-related protection systems to monitor coolant
flow, over temperature, etc. The main beam parameters of concern are density
(sometimes derived from a beam profile measurement), beam size, beam position,
and beam current (from which beam power can be calculated). Beam-related
target protection involves control over these parameters, monitoring these
parameters, and turning off the beam if these parameters move outside
allowable limits. Types of beam-related target protection include:

  * •

rapid and automatic beam turn off for off-normal beams;

  * •

locking down equipment and operating control parameters to prevent accidental
changes;

  * •

alarms that alert operators when certain parameters move outside of pre-
defined limits;

  * •

collimators to partially intercept off-normal beams;

  * •

beam transport designs that avoid high sensitivity—e.g. do not want to live on
the edge, do not want to be very sensitive to small magnet changes;

  * •

target designs that can tolerate off-normal beam for a short period of time
(e.g. most SNS equipment is designed such that it can withstand two full-power
off-normal beam pulses).

Figure 14:  Some typical turn-off times for the SNS case

###  0.7.1 Protection by rapid beam turn off

The time it takes to automatically turn off the beam varies depending on the
method used. Figure 14 shows some typical turn-off times for the SNS case.

The fastest method to turn off the beam is through the MPS. After all, that is
the purpose of the MPS. A typical time to turn off the beam at SNS is 20 to 30
$\mu$s for trips that originate in the linac, for example if a beam loss
monitor exceeds a pre-determined threshold. Trips that originate further
downstream, in the ring or in the ring-to-target beam transport, can take a
few $\mu$s longer. An interesting fine point is that if the trip occurs
upstream of or within the accumulator ring, the MPS can stop beam injection,
but the beam that has already been injected into the ring will be sent to the
target. The beam intensity can be anywhere from small to full intensity. If
the trip occurs downstream of the ring, it can only prevent the next pulse
from beginning the accumulation process.

For example, at SNS, PLCs monitor magnet current transformers on the power
supplies for the last few magnets in the transport line to the target. If the
magnet current strays outside pre-determined limits, the beam will be
automatically turned off. The time for a PLC trip to turn off the beam is
about 50 ms. As an additional precaution, PLC trips also turn off the high
voltage to the ion source and turn off the high voltage to the first few
klystrons.

The third-fastest method is faults detected by the control computers. The
computers interface to the MPS to turn off the beam. For example, at SNS,
computers monitor all the magnet power supply currents in the beam transport
to the target (as read by the power supply controller—not by separate current
transformers). The estimated time to turn off the beam is about 1 s. For
example, Fig. 15 shows a partial list of these magnets and their high- and
low-trip limits.

Figure 15:  Some of the power supplies in the SNS beam transport from the ring
to the target, and their trip limits that are monitored by the control
computer.

The fourth-fastest method is the alarm system, which alerts the operators to
off-normal conditions. The alarms are generated by the control system, and
require operator action to turn off the beam or to correct the off-normal
condition. A strength of this method is that the control system can monitor a
huge number of parameters, and also values derived from those parameters. For
example, at SNS, the control system alarms if the beam position on target
strays outside of pre-determined limits. The estimated time to correct the
alarm state is 10 to 30 s, or possibly longer, depending on the operator
reaction time and the complexity of the alarm.

Beam loss monitors (BLMs) are a standard part of machine protection for any
beam transport system, and their purpose is primarily to protect beam line
components. But BLMs can also be used to protect a target. The idea is that a
quadrupole magnet change that can cause the beam to be too small or too large
at the target will sometimes also cause beam loss upstream of the target,
because the mis-focused beam could be large enough somewhere along the beam
line that the beam tails will strike the beam pipe walls. BLMs can also
protect against deviations from the nominal beam trajectory that could result
in a bad position on the target, if the deviation is large enough to cause
beam loss. At SNS, BLMs are located, and thresholds are set, such that some of
the quadrupole magnet changes that can result in off-normal beams at the
target will trip the BLMs. Figure 16 shows, for each quadrupole magnet in the
beam line from the ring to the target, how much the beam size on the target
changes for a 5% change in the quadrupole magnet gradient. It also shows the
greatest amount of change anywhere in the beam line. Figure 17 shows the six
of 19 quadruple magnets that will cause a BLM trip (beam loss greater than
0.1%) before causing dangerous beam parameters at the target.

Figure 16:  These plots show, for each quadrupole magnet in the beam line from
the ring to the target, how much the beam size on the target changes for a 5%
change in the quadrupole magnet gradient (blue line). Also shown is the
maximum beam size fractional change along the beam line (red line). Figure
reproduced from Ref. [20]. Figure 17:  The six of 19 quadrupole magnets that
will cause a BLM trip (beam loss greater than 0.1%) before causing dangerous
beam parameters at the target. Figure reproduced from Ref. [20].

###  0.7.2 Protection by equipment lock down

Another method to protect the target is to lock down certain control
parameters. The idea is to administratively control critical hardware set
points to prevent inadvertent changes. An example of this method is that at
SNS the operators set a gate generator at the ion source that limits the
possible beam pulse lengths, which limits the beam power on the target.

###  0.7.3 Protection by collimation

High beam power facilities often have collimators in the last part of the beam
transport leading to the target, to protect against large beam position
variations and overly large beam sizes, and to ensure that the beam hits the
central region of the target. For example, SNS has a collimator immediately
upstream of the target. The $27.9\times 12.7$ cm2 aperture is smaller than the
face of the target, but slightly larger than the nominal beam size (90% of
beam must fit within a $20\times 7$ cm2 rectangle). SNS also has collimators
in the beam transport line that can intercept a portion of the beam that does
not receive the nominal kick angle from the extraction kickers, thus
protecting the target from some of the mis-kicked beam that would consequently
arrive off centre at the target. Figure 18 shows the collimator that is
embedded in the bulk shielding surrounding the SNS target, Fig. 19 shows the
positions of the beam line collimators, and Fig. 20 shows simulations of the
beam trajectory deviations due to extraction kickers that do not fire
properly.

Figure 18:  The collimator and shielding just upstream of the SNS neutron
production target. The PBW label refers to the Proton Beam Window. Figure 19:
The downstream portion of the SNS RTBT beam line. The blue arrows show the
positions of the RTBT collimators. Figure 20:  Simulations of beam trajectory
deviations caused by missed kicks in the ring extraction kicker system. The
thin lines indicate trajectory deviations for individual missed kicks, and the
thick lines indicate the maximum and minimum trajectory deviations for the
case of two simultaneous missed kicks. Figure reproduced from Ref. [21].

###  0.7.4 Protection by interlocking on off-normal measured beam parameters

This protection category involves continuously monitoring critical beam
parameters, and then automatically turning off the beam if the parameters
exceed pre-determined thresholds. Table 2 shows some of the more common beam
instrumentation used for this purpose.

Table 2: Some example target interlocks Interlock |  Some accelerators that use this as a target interlock
---|---
Beam loss monitors |  Everybody
Beam current monitors |  Everybody
Beam position monitors |  PSI and FNAL have automatic beam centring based on beam position monitor signals
Harp profile monitors |  ISISa, LANSCEa, J-PARCa, SNSa
Ionization profile monitors |  PSI
Glowing screen just upstream of target (profile monitor) |  PSI
Halo thermocouples |  SNS
Amount of beam intercepted by collimator |  PSI
a No interlock at present but this is under consideration

###  0.7.5 Protection by beam line design

It is usually desirable to design the beam transport line to avoid high
sensitivity to beam parameters at the target. For example, small adjustments
to dipole magnets should not cause large position changes on the target, and
small adjustments to quadrupole magnets should not cause large beam size
changes on the target. At SNS, one way this design practice was implemented
was to require the phase advance from the extraction kickers to the target to
be approximately a multiple of $\pi$. In this way, a missed kick from an
extraction kicker will not cause a large position change on the target. This
effect is illustrated in Fig. 20.

##  0.8 Some example target protection implementations

###  0.8.1 J-PARC target protection

The J-PARC design parameters call for a 1 MW, 3 GeV, 25 Hz proton beam to be
delivered to a liquid mercury target. To protect the target [22], there is a
waveform monitor on each of the eight extraction kickers. A PLC monitors some
of the quadrupole and dipole magnet power supply currents. A PLC also monitors
the beam profile measurement device upstream of the target, and it is possible
to connect this to the machine protection system, but presently this is not
the case. A recent addition to the MPS interlock system is a fast beam current
monitor to shut off the beam if the beam current exceeds a pre-defined limit.

###  0.8.2 FNAL NuMI target protection

The NuMI neutrino production target at FNAL now operates with 120 GeV,
approximately 400 kW protons, <1 Hz, with an upgrade to 700 kW planned in
2015. A later upgrade to >1 MW is planned. Prior to beam extraction from the
main injector (MI), more than 250 different inputs to the beam permit system
are checked [23]. The beam is sent to the abort dump if inputs are not
correct. The inputs checked include the following: beam position and angle for
the MI extraction channel, possible excessive residual MI beam in the kicker
beam gap, the extraction kicker status, proper NuMI power supply flat-top
values, the beam readiness of the target station and absorbers, all beam loss
monitor readings from the previous extraction, and the previous pulse position
and trajectory at the target. The beam delivery system also has an automatic
beam steering system to keep the beam centred on the target. If a beam pulse
is >1.5 mm off centre, the beam is automatically turned off.

###  0.8.3 ISIS target protection

The ISIS facility at RAL operates two targets, TS1 and TS2. TS1 operates with
800 MeV, 160 kW, 40 Hz protons, and TS2 operates with 800 MeV, 40 kW, 10 Hz
protons. A beam halo monitor (comprising eight equally spaced thermocouples in
the penumbra of the beam $\sim$100 mm upstream of the target) indicates mis-
focused and mis-steered beams. If the thermocouple signals exceed pre-
determined trip thresholds, the protection system turns off the beam. The
protection system [24] is capable of turning off the machine within 2 ms and
hence inhibits the next beam pulse from the ion source in the 50 Hz
synchrotron cycle. There is also a steering servo system to accommodate
variations in the upstream beam. The servo and trip system is slow, running at
$\sim$2 s. On TS2 there is also a harp profile monitor which sits permanently
in the beam, but this is not interlocked. A similar harp is planned for TS1 as
part of an upgrade project.

##  0.9 SNS target tune up

To provide a detailed example of target protection, we will now discuss the
SNS case. We will be referring to the beam instrumentation shown in Fig. 21.
The initial beam tuning to the target is performed at low intensity, low
power, 1 Hz (<1 kW) beam. Once the beam position has been adjusted to centre
the beam on the target, the beam power is slowly increased to full intensity,
but still at 1 Hz ($\sim$23 kW). The beam intensity is now high enough that an
image can be seen on the target imaging system (TIS), which is based on a
light-emitting coating on the surface of the target, and hence gives a direct
measurement of beam position and beam distribution at the target. The TIS is
used to check the beam centring, and then the beam size is measured not with
the TIS (due to unresolved discrepancies with the measurement technique about
to be discussed), but with the four wire scanners and one harp located at
various positions in the ring to target beam transport (RTBT) beam line, as
shown in Fig. 21. The beam size at the target is determined by fitting the
measured rms beam sizes with an online model, and varying the simulated beam
parameters until the best fit to the measured rms beam sizes is obtained. The
model is then used to extrapolate the beam size at the target. An example fit
is shown in Fig. 22, and Fig. 23 shows the corresponding output table with the
extrapolated beam size at the target.

The same online model is then used to extrapolate the peak beam density,
measured with the harp, to the target. The critical beam parameters for the
target have now been measured. To meet requirements:

  * •

the beam must be centred on the target to within 6 mm horizontal and 4 mm
vertical;

  * •

the rms beam size must be <49 mm horizontal and <17 mm vertical;

  * •

the peak density on target must be less than $2\times 10^{16}$ protons/m2;

  * •

the peak density on the proton beam window must be less than $2.9\times
10^{16}$ protons/m2.

Once the measured beam parameters have been determined to meet requirements,
the beam delivery system is locked down. This is required before increasing
the beam power above 100 kW. The lock down is accomplished by engaging the
following interlock thresholds.

  * •

The last five quadrupole magnet power supplies are monitored by a PLC with
current limits set to $\pm$7%.

  * •

Both large dipoles are monitored by a PLC with current limits set to $\pm$2 A.

  * •

The last two horizontal dipole correctors and last two vertical dipole
correctors are monitored by a PLC with current limits set to $\pm$5 A.

  * •

The injection kicker waveform monitor system is engaged. If waveforms stray
outside of pre-defined windows, the MPS trips the beam.

  * •

The extraction kicker waveform monitor system is engaged. If waveforms stray
outside of pre-defined windows, the MPS trips the beam.

  * •

All RTBT magnets are monitored by the control computer with current limits set
to $\pm$5% on quadrupole magnets, $\pm$5% on large dipole magnets, and
$\pm$0.5 A on dipole corrector magnets.

  * •

The beam power limit is set into the control system, and monitored by the
computer.

The beam pulse length and beam duty factor are also locked down, but in a way
that does not cause a beam trip. It simply prevents the timing system hardware
from being set in a way that could exceed a pre-determined beam power. A
screen shot of this system is shown in Fig. 24.

Figure 25 shows a screen shot of the magnet currents that are monitored by a
PLC. The PLC monitors those magnets that have the biggest impact on beam size
and position at the target. The control system computers also monitor these
magnets, as well as many more, as shown in Fig. 26, but the PLC is a more
robust and reliable system compared to the control computer.

There are four horizontal and four vertical injection kickers to paint the
beam distribution in the ring. An example set of waveforms is shown in Fig.
27. If for some reason the waveforms were accidentally set to paint a smaller
beam, it would result in an excessively high beam density on the target. Or,
the kicker system may fail in a way that could also cause the beam size to be
too small. To protect against these possibilities an injection kicker waveform
monitor system, based on commercially available oscilloscopes, watches the
readback waveforms and automatically trips off the beam if the waveforms stray
outside pre-determined envelopes.

Some interlocks are always in effect, not just when the beam power is greater
than 100 kW. Examples include the beam loss monitors (direct MPS interlock),
the proton beam window halo thermocouples (shown in Fig. 5, and interlocked
through a PLC connected to the MPS), and the target protection system, which
monitors parameters such as mercury flow, water cooling, etc (PLC interlock).

The control system computers monitor the following parameters and alert the
operators if they stray outside of pre-set bounds:

  * •

beam position on target calculated by the target imaging system;

  * •

the beam density at the harp 9.52 m upstream of the target. The beam density
at the target should be proportional to the beam density at the harp;

  * •

the beam size at the harp 9.52 m upstream of the target. The beam size at the
target should be proportional to the beam size at the harp;

  * •

beam power;

  * •

beam centring estimated from the proton beam window halo thermocouples
(top–bottom, left–right);

  * •

some magnet currents.

The last three items on the list are also interlocked by other methods as
discussed above, but in this case the thresholds are more loosely set to give
the operators a chance to correct the problem before it becomes severe enough
to trip off the beam. The alarm summary screen is shown in Fig. 28.

Figure 21:  A drawing of the SNS beam line upstream of the target, showing the
locations of the profile and position monitors. Image reproduced from Ref.
[25]. Figure 22:  Screen shot of the application that fits a model to measured
beam sizes to extrapolate the beam size at the target. The points show the
measured beam sizes and the lines show the model beam sizes. Figure 23:
Another screen shot of the application that extrapolates the beam size on the
target. The resultant size predictions at the beam window and target are
circled. Figure 24:  A screen shot of the application showing the DG535 pulse
generator settings that lock down the maximum possible beam pulse width.
Figure 25:  A screen shot showing the PLC limits on the last magnets upstream
of the target Figure 26:  A screen shot showing the magnets monitored by the
control system to lock down the beam transport to the target. Figure 27:
(Colour) A screen shot showing the injection kicker waveforms. The region
marked by the rectangle and arrow shows the portion of the waveform that is
monitored for one of the kickers (note that all eight are monitored). Figure
28:  A screen shot of the alarm summary screen and example screen shots of the
beam profile measured by the harp and the beam distribution measured by the
TIS. The arrows show how the measured beam parameters are monitored by the
alarm system.

##  0.10 Summary

In summary, high-power targetry is an active field that is continuously
advancing. The highest power targets require tight control over the beam
position, density, and distribution, and even then the target lifetimes can be
short (e.g. 6 months at SNS). Machine protection systems must monitor these
beam parameters and quickly activate interlocks if they stray outside of
limits. Protection can include monitoring equipment set points, equipment
readback parameters, beam parameters, locking down certain controls, beam line
design, and collimator systems.

## Acknowledgements

ORNL is managed by UT-Battelle, LLC, under contract DE-AC05-00OR22725 for the
U.S. Department of Energy.

Notice: This manuscript has been authored by UT-Battelle, LLC, under Contract
No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains, and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide licence to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes.

## References

  * [1] J. Haines, U.S. Particle Accelerator School, High Power Beam Targets Class, Vanderbilt University, Nashville, Tennessee, 20 January 2009.
  * [2] B. Reimer, "Spallation Source Facilities," 5th High Power Targetry Workshop, Fermilab, 20–23 May 2014. https://indico.fnal.gov/conferenceDisplay.py?confId=7870.
  * [3] D. Haynes, "Introduction and Overview of the ISIS Target Systems," ISIS and SNS Bilateral Workshop, Oak Ridge, Tennessee, 25–26 June 2013.
  * [4] T. McManamy, U.S. Particle Accelerator School, High Power Beam Targets Class, Vanderbilt University, Nashville, Tennessee, 20 January 2009.
  * [5] R. Werbeck, "LANSCE Short-Pulse Target Operation," Workshop on High-power Targetry for Future Accelerators, Long Island, New York, 8–12 September 2003.
  * [6] W. Wagner et al., "The SINQ solid spallation target ? operational experience and recent improvements ," High Power Targetry Workshop, Malmo, Sweden, 2–6 May 2011. Also http://www.hep.princeton.edu/mumu/target/Wagner/Wagner_050311.pdf.
  * [7] N. Mokhov, "Beam-Material Interactions," Joint US–CERN Accelerator School on Machine Protection, Newport Beach, CA, USA, 5–14 November 2014.
  * [8] J. Hylen, "Survey of Target Facility Landscape: Neutrino Beam Facilities," 5th High Power Targetry Workshop, Fermilab, 20–23 May 2014.
  * [9] P. Bricault, "Radioactive Ion Beam Facilities ? High Power Target, Current status and Future Directions," 5th High Power Targetry Workshop, Fermilab, 20–23 May 2014. Also https://indico.fnal.gov/getFile.py/access?contribId=91&sessionId=0&resId=0&materialId=slides&
confId=7870.

  * [10] P. Hurh, Third High Power Target Workshop, Bad Zurach, Switzerland, September 2007.
  * [11] T. McManamy, "Operational Experience with the SNS Target Systems and Upgrade Plans," Workshop on Applications of High Intensity Proton Accelerators, Fermi National Accelerator Laboratory, Batavia, IL, USA, 19–21 October 2009.
  * [12] S. Meigo, "Development of beam flattening system using non-linear beam optics at J-PARC," Eleventh International Topical Meeting on Nuclear Applications of Accelerators, Bruges, Belgium, 5–8 August 2014.
  * [13] S. Meigo et al., Beam flattening system based on non-linear optics for high power spallation neutron target at J-PARC, Proc. IPAC2014, Dresden, Germany, 15–20 June 2014.
  * [14] H.D. Thomsen, A.I.S. Holm, and S.P. Møller, "A linear beam raster system for the European Spallation Source," Proc. IPAC2013, Shanghai, China, May 2013, p. 70.
  * [15] http://j-parc.jp/en/topics/20130812Accident_Report.html; http://j-parc.jp/en/topics/HDAccident
20131217.pdf

  * [16] L. Jones, PASI Working Group Report by D.M. Jenkins, April 2013; http://pasi.org.uk/images/2/2c/DavidJenkins_PASI_meeting_28_Feb_13.pdf.
  * [17] K. Thomsen and P.A. Schmelzbach, "A dedicated beam interrupt system for the safe operation of the Megapie Liquid Metal Target," Utilisation and Reliability of High Power Proton Accelerators (HPPA5), Workshop Proc., Mol, Belgium, 6–9 May 2007.
  * [18] K. Thomsen, "VIMOS, near-target beam diagnostics for MEGAPIE," NIM A 575(3) (2007) 347. http://dx.doi.org/10.1016/j.nima.2007.03.011
  * [19] W. F. Sommer et al., "Failure analysis of a radio-activated accelerator component," Practical Failure Anal. 3(1) (2003) 71. http://dx.doi.org/10.1007/BF02717412
  * [20] J. Holmes, "Quadrupole Strength Limits in the RTBT," SNS Tech. Note 162, 2005.
  * [21] N. Catalan-Lasheras and D. Raparia, "The Collimation System of the SNS Transfer Lines," Proc. of the 2001 Particle Accelerator Conference, Chicago, p. 3263. http://accelconf.web.cern.ch/AccelConf/p01/INDEX.HTM
  * [22] M. Kinsho, Private communication, 2014.
  * [23] S. Childress, "NUMI Proton Beam Diagnostics and Control: Achieving 2 Megawatt Capability," Proc. HB2008 Workshop, Nashville, TN, USA, 25–29 August 2008, p. 475.
  * [24] D. Adams, Private communication, 2014.
  * [25] M. Plum, "SNS Injection and Extraction Systems - Issues And Solutions," Proc. HB2008 Workshop, Nashville, TN, USA, 25–29 August 2008, p. 268.


# On the dynamics and integrability
of the Ziegler pendulum

Ivan Yu. Polekhin Electronic address<EMAIL_ADDRESS>Steklov
Mathematical Institute of Russian Academy of Sciences,
Moscow, Russia
<EMAIL_ADDRESS>
## Abstract

We prove that the Ziegler pendulum — a double pendulum with a follower force —
can be integrable, provided that the stiffness of the elastic spring located
in the pivot point of the pendulum is zero. We show that the integrability of
the system follows from the existence of two-parameter families of periodic
solutions. We explain the mechanism for the transition from the integrable to
the chaotic dynamics and show that this mechanism is different from the
typical phenomena observed in conservative Hamiltonian systems. The case when
the stiffnesses of the both springs are non-zero is studied numerically. We
show that the regular dynamics co-exists with the chaotic dynamics.

Key words:  first integrals, non-potential force, system on a torus, non-
integrability, periodic solutions.

##  1 Introduction

The Ziegler pendulum is a well known non-conservative mechanical system with
many real-life applications. This system is mostly known for the ‘paradox’ of
loss of stability in the presence of friction, which was discovered 70 years
ago [1] (see also [2, 3]). More detailed exposition of this and related
systems and the history of the problem can be found in [4]. The full
bibliography of works on the Zielger pendulum and other non-conservative
systems is too extensive to be presented here and we also refer the reader to
book [4]. However, we would like to note that, for the best of our knowledge,
the dynamics in the large and integrability of the Zielger pendulum has not
been studied before in the literature. At the same time, there are many papers
on the dynamics of the classical double pendulum, where it is shown
experimentally [5], numerically [6], and analytically [7] by means of the
Poincare-Melnikov method that this system is chaotic. First non-trivial
results on the non-integrability of systems in non-potential force fields have
been obtained in [8]. The problem of integrability of the Ziegler pendulum has
been also formulated in [8].

In the paper we study the dynamics of a planar double pendulum in the presence
of a follower force of constant magnitude. We assume that this follower force
is acting along the rod that is not connected to the fixed pivot point. There
is no gravity force in the system, yet there are two springs of linear
stiffness located in the joints of the double pendulum.

The paper contains both analytical and numerical results. By means of an
analytical approach we study the case when the stiffness of the spring located
in the pivot point is zero. Thanks to the rotational symmetry, we can consider
the reduced system of the third order. We prove that in the reduced system
there exists a two-parametric family of periodic solutions. It means that
locally there exist two first integrals of the system and one can consider
this system as integrable, at least for some initial conditions. We explain
the mechanism of loss of periodicity of solutions associated with the change
of parameters or initial conditions. We show how the birth of non-periodic
solutions leads to the non-integrability of the system.

The case when the stiffnesses of the springs are non-zero is studied
numerically. We again obtain the co-existence of two types of solutions in the
phase space: regular solutions which lie on two-dimensional surfaces and
chaotic trajectories. We show this by means of the Poincare sections and
Lyapunov characteristic exponents.

##  2 Equations of motion

Let us consider a planar double mathematical pendulum. One of its absolutely
rigid massless rod of length $l_{2}$ can rotate freely in the plane about a
fixed point $O$. Two other rigid massless rods of lengths $l_{1}$ and $l_{3}$
are connected to the free end $A$ of the first rod $OA$. Let us denote the
ends of these rods by $B$ and $C$, correspondingly (Fig. 1). We assume that
points $A$, $B$, and $C$ always remain on a straight line when the pendulum
moves. There is a point mass $m_{2}$ in point $A$ and point masses $m_{1}$ and
$m_{3}$ in points $B$ and $C$, correspondingly. There is a force $F$ of
constant magnitude acting along segment $BC$. Two springs of stiffnesses
$k_{1}$ and $k_{2}$ are located in points $A$ and $O$, correspondingly.

Figure 1: The Ziegler pendulum.

Let $\varphi_{2}$ be a generalized coordinate, the angle between a constant
direction in the plane and segment $OA$ and $\varphi_{1}$ be the angle between
$OA$ and $BC$. The kinetic energy of the system can be expressed in these
coordinates as follows:

$\displaystyle\begin{split}T=&\frac{m_{1}}{2}\left(\dot{\varphi}_{2}^{2}l_{2}^{2}+(\dot{\varphi}_{1}+\dot{\varphi}_{2})^{2}l_{1}^{2}+2l_{1}l_{2}\dot{\varphi}_{2}(\dot{\varphi}_{1}+\dot{\varphi}_{2})\cos\varphi_{1}\right)+\frac{m_{2}}{2}\dot{\varphi}_{2}^{2}l_{2}^{2}+\\\
&\frac{m_{3}}{2}\left(\dot{\varphi}_{2}^{2}l_{2}^{2}+(\dot{\varphi}_{1}+\dot{\varphi}_{2})^{2}l_{3}^{2}-2l_{2}l_{3}\dot{\varphi}_{2}(\dot{\varphi}_{1}+\dot{\varphi}_{2})\cos\varphi_{1}\right).\end{split}$
(1)

The potential energy:

$\displaystyle\begin{split}\Pi=\frac{k_{1}\varphi_{1}^{2}}{2}+\frac{k_{2}\varphi_{2}^{2}}{2}.\end{split}$
(2)

The components of generalized force:

$\displaystyle\begin{split}Q_{1}=0,\quad
Q_{2}=-Fl_{2}\sin\varphi_{1}.\end{split}$ (3)

From (1), (2), and (3) we obtain the Lagrange equations, which can be
presented as follows

$\displaystyle\begin{split}&A_{11}\ddot{\varphi}_{1}+A_{12}\ddot{\varphi}_{2}=r_{1},\\\
&A_{21}\ddot{\varphi}_{1}+A_{22}\ddot{\varphi}_{2}=r_{2},\end{split}$ (4)

where

$\displaystyle\begin{split}&A_{11}=m_{1}l_{1}^{2}+m_{3}l_{3}^{2},\\\
&A_{12}=A_{21}=m_{1}l_{1}^{2}+m_{3}l_{3}^{2}+m_{1}l_{1}l_{2}\cos\varphi_{1}-m_{3}l_{2}l_{3}\cos\varphi_{1},\\\
&A_{22}=m_{1}l_{2}^{2}+m_{1}l_{1}^{2}+m_{2}l_{2}^{2}+m_{3}l_{2}^{2}+m_{3}l_{3}^{2}+2m_{1}l_{1}l_{2}\cos\varphi_{1}-2m_{3}l_{2}l_{3}\cos\varphi_{1},\\\
&r_{1}=-k_{1}\varphi_{1}-m_{1}l_{1}l_{2}\dot{\varphi}_{2}^{2}\sin\varphi_{1}+m_{3}l_{2}l_{3}\dot{\varphi}_{2}^{2}\sin\varphi_{1},\\\
&r_{2}=-Fl_{1}\sin\varphi_{1}-k_{2}\varphi_{2}+m_{1}l_{1}l_{2}\dot{\varphi}_{1}(\dot{\varphi}_{1}+2\dot{\varphi}_{2})\sin\varphi_{1}-m_{3}l_{2}l_{3}\dot{\varphi}_{1}(\dot{\varphi}_{1}+2\dot{\varphi}_{2})\sin\varphi_{1}.\end{split}$

If $k_{2}=0$ then the right-hand side does not depend on variable
$\varphi_{2}$ and the system of the third order for variables $\varphi_{1}$,
$\dot{\varphi}_{1}=v_{1}$, and $\dot{\varphi}_{2}=v_{2}$ can be considered
independently

$\displaystyle\begin{split}&\dot{\varphi}_{1}=v_{1},\\\
&\dot{v}_{1}=\left(r_{1}-\frac{A_{12}}{A_{22}}r_{2}\right)\left(A_{11}-\frac{A_{12}A_{21}}{A_{22}}\right)^{-1},\\\
&\dot{v}_{2}=\left(r_{2}-\frac{A_{21}}{A_{11}}r_{1}\right)\left(A_{22}-\frac{A_{12}A_{21}}{A_{11}}\right)^{-1}.\end{split}$
(5)

An important property of system (4) is the existence of a smooth invariant
measure. Indeed, it is a system with circulatory forces, i.e. the Lagrange
equations have the form

$\frac{d}{dt}\frac{\partial L}{\partial\dot{q}_{i}}-\frac{\partial L}{\partial
q_{i}}=Q_{i}(q),\quad i=1,2.$

Since the components of the generalized force depend only on the generalized
coordinates and do not depend on the generalized momenta, then in the
Hamiltonian coordinates this system has an invariant measure of constant
density.

Let us briefly discuss various definitions of integrability which can be
applied to systems (4) and (5). First, if system (4) is Hamiltonian (for
instance, for $F=0$), then it is natural to call the system integrable if the
Liouville-Arnold theorem holds. Second, if this system is not Hamiltonian, we
can consider integrability in the sense of the Jacobi theorem on the last
multiplier (see, for instance, [9]). We say that system (4) is integrable if
there exist two independent first integrals. Similarly, we say that system (5)
is integrable if there exist two first integrals for this system. In this case
the integrability follows from the classical theorem on the integrability for
general ODE systems: a system in $\mathbb{R}^{n}$ is integrable if it has
$n-1$ first integrals.

##  3 Cases $F=0$ and $m_{1}l_{1}=m_{3}l_{3}$

Let us consider two cases in which the dynamics of system (5) (for $k_{2}=0$)
can be fully understood. First, if there is no follower force, then the system
is Hamiltonian and there exist two first integrals: the energy integral
$H=T+\Pi$ and the cyclic integral (the integral of angular momentum)

$\displaystyle\begin{split}K=\frac{\partial
T}{\partial\dot{\varphi}_{2}}=&\frac{m_{1}}{2}\left(2\dot{\varphi}_{2}l_{2}^{2}+2(\dot{\varphi}_{1}+\dot{\varphi}_{2})l_{1}^{2}+2l_{1}l_{2}\cos\varphi_{1}(\dot{\varphi}_{1}+2\dot{\varphi}_{2})\right)+m_{2}\dot{\varphi}_{2}l_{2}^{2}+\\\
&\frac{m_{3}}{2}\left(2\dot{\varphi}_{2}l_{2}^{2}+2(\dot{\varphi}_{1}+\dot{\varphi}_{2})l_{3}^{2}-2l_{2}l_{3}(\dot{\varphi}_{1}+2\dot{\varphi}_{2})\cos\varphi_{1}\right).\end{split}$
(6)

Therefore, the system is integrable in the sense of the Liouville-Arnold
theorem, as a general system of ordinary differential equations (system (5)
defined in a three-dimensional space has two first integrals), and in the
sense of the Jacobi theorem on the last multiplier. The latter follows from
the fact that the original system (4) has an invariant measure for any $F$.
The density of this measure is constant if we rewrite our system in the
Hamiltonian coordinates.

Another case in which system (5) is simplified and can be integrated in
quadratures is when $m_{1}l_{1}=m_{3}l_{3}$. This equality allows one to
consider two separate independent systems:

$\displaystyle\begin{split}&\ddot{\varphi}_{1}=\frac{Fl_{1}\sin\varphi_{1}-k_{1}\varphi_{1}}{(m_{1}+m_{2}+m_{3})l_{2}^{2}},\\\
&v_{2}(t)=v_{2}(0)-\int\limits_{0}^{t}\ddot{\varphi}_{1}+\frac{k_{1}}{m_{1}l_{1}^{2}+m_{3}l_{3}^{2}}\varphi_{1}\,dt.\end{split}$
(7)

This system is integrable by the Jacobi theorem: the first integral is the
energy integral of the first equation, the phase volume is an invariant
measure, if we present the system in the Hamiltonian form.

The periodic solutions that intersect the plane $\varphi_{1}=0$ in two points,
are of particular interest for our considerations. We show that such solutions
exist in the both considered cases (Fig. 2).

(a) $F=0$

(b) $m_{1}l_{1}=m_{3}l_{3}$

Figure 2: Examples of families of periodic solutions intersecting plane
$\varphi_{1}=0$ in two points.

######  Proposition 1.

For $F=0$ and large $k_{1}$ there exists a two-parameter family of periodic
solutions of system (5) which intersect the plane $\varphi_{1}=0$
transversely.

###### Proof.

Find an expression for $\dot{\varphi}_{2}$ from (6):

$\dot{\varphi}_{2}=\frac{K-V(\varphi_{1})\dot{\varphi}_{1}}{U(\varphi_{1})},$

where

$\displaystyle\begin{split}&V(\varphi_{1})=m_{1}l_{1}^{2}+m_{3}l_{3}^{2}+m_{1}l_{1}l_{2}\cos\varphi_{1}-m_{3}l_{2}l_{3}\cos\varphi_{1},\\\
&U(\varphi_{1})=m_{1}l_{1}^{2}+m_{1}l_{2}^{2}+m_{2}l_{2}^{2}+m_{3}l_{2}^{2}+2m_{1}l_{1}l_{2}\cos\varphi_{1}-2m_{3}l_{2}l_{3}\cos\varphi_{1}.\end{split}$

Substitute this expression for $\dot{\varphi}_{2}$ into the energy integral

$\displaystyle\begin{split}H(\varphi_{1},\dot{\varphi}_{1})=&\frac{1}{2}l_{2}^{2}(m_{1}+m_{2}+m_{3})\left(\frac{K-V\dot{\varphi}_{1}}{U}\right)^{2}+\frac{1}{2}(m_{1}l_{1}^{2}+m_{3}l_{3}^{2})\left(\frac{K+(U-V)\dot{\varphi}_{1}}{U}\right)^{2}+\\\
&\left(\frac{K-V\dot{\varphi}_{1}}{U}\right)\left(\frac{K+(U-V)\dot{\varphi}_{1}}{U}\right)\cos\varphi_{1}(m_{1}l_{1}l_{2}-m_{3}l_{2}l_{3})+\frac{k_{1}\varphi_{1}^{2}}{2}.\end{split}$

Let us prove that there is an extremum on the axis $\varphi_{1}=0$. Since
$U^{\prime}(0)=V^{\prime}(0)=0$, then for $\varphi_{1}=0$, we have

$\frac{\partial H}{\partial\varphi_{1}}=0.$

From eqution

$\frac{\partial H}{\partial\dot{\varphi}_{1}}=0$

we obtain

$\dot{\varphi}_{1}=\frac{VKl_{2}^{2}(m_{1}+m_{2}+m_{3})+K(V-U)(m_{1}l_{1}^{2}+m_{3}l_{3}^{2})+(2V-U)K(m_{1}l_{1}l_{2}-m_{3}l_{2}l_{3})}{V^{2}l_{2}^{2}(m_{1}+m_{2}+m_{3})+(U-V)^{2}(m_{1}l_{1}^{2}+m_{3}l_{3}^{2})-2V(U-V)(m_{1}l_{1}l_{2}-m_{3}l_{2}l_{3})}.$

It can be shown that the denominator of the fraction is strictly positive for
all values of parameters. In particular, we obtain that the first minor of the
Hessian matrix is positive. Non-diagonal elements of this matrix equal zero
and the second diagonal element is positive provided that $k_{1}$ is
sufficiently large.

Finally, we have that for all values of the energy integral which is close to
the critical value of function $H$ at the point of local minimum, the
trajectory is periodic and intersects the plane $\varphi_{1}=0$ transversely.
Moreover, this result holds for all values of first integral $K$.

∎

######  Proposition 2.

For $k_{1}>0$ system (7) has a two-parameter family of periodic solutions
which intersect the plane $\varphi_{1}=0$ transversely.

###### Proof.

The first equation of system (7) is the Lagrange equation with the following
potential energy

$\Pi=\frac{Fl_{1}\cos\varphi_{1}+\frac{1}{2}k_{1}\varphi_{1}^{2}}{(m_{1}+m_{2}+m_{3})l_{2}^{2}}.$

If the value of total energy is sufficiently large, then the phase
trajectories encircle the point of origin $\varphi_{1}=0$,
$\dot{\varphi}_{1}=0$ and these trajectories are symmetric with respect to the
coordinate axes. Therefore, $\varphi_{1}(t)$ is a periodic function and its
mean value equals zero. We obtain that $v_{2}(t)$ is also a periodic function
with the same period. ∎

######  Proposition 3.

For $k_{1}=0$ and $F<0$ system (7) a two-parameter family of periodic
solutions which intersect the plane $\varphi_{1}=0$ transversely.

###### Proof.

For the first equation of the system, point $\varphi_{1}=0$,
$\dot{\varphi}_{1}=0$ is a local minimum. ∎

Note that the Zielger pendulum can be considered as a Hamiltonian system for
some values of the parameters [10]. However, for the case $k_{2}=0$ it follows
from the conditions presented in [10] that $m_{1}l_{1}=m_{3}l_{3}$ and
$k_{1}=0$. Therefore, from (4) we obtain a trivial system, which is less
general than (7).

##  4 Integrability

For the case when $k_{2}=0$ it is natural to call the original system
integrable if the reduced system (5) has two first integrals. In the general
case, for which it is also allowed that $k_{2}\neq 0$, it is also natural to
call the system with two first integrals integrable. For any $k_{2}$ we have a
system defined on $\mathbb{T}^{2}\times\mathbb{R}^{2}$ and this system has an
invariant measure. Therefore, this system is integrable if it has two first
integrals, similarly to Hamiltonian systems with two degrees of freedom and to
the case $k_{2}=0$.

Let us prove that system (5) for $k_{2}=0$ remains integrable for $|F|$ or
$|m_{1}l_{1}-m_{3}l_{3}|$ sufficiently small. The key step in the proof is
given by the following lemma.

######  Lemma 1.

Let functions $\varphi_{1}(t)$, $v_{1}(t)$, $v_{2}(t)$ be a solution of system
(5). Then functions $-\varphi_{1}(-t)$, $v_{1}(-t)$, $v_{2}(-t)$ are also a
solution.

###### Proof.

The proof follows from the substitution of given functions into the system.
Terms $A_{11}$, $A_{12}$, $A_{22}$ do not change after substitution
$\varphi_{1}\mapsto-\varphi_{1}$. Terms $r_{1}$ and $r_{2}$ change their
signs. At the same time, the left-hand side also changes its sign. ∎

From the lemma we obtain the following

######  Proposition 4.

Let $\varphi_{1}(t)$, $v_{1}(t)$, $v_{2}(t)$ be a solution which intersects
the plane $\varphi_{1}=0$ in two points. Then this solution is periodic.

###### Proof.

A periodic solution can be obtained from the segment between the two points
using the symmetry with respect to the plane $\varphi_{1}=0$. ∎

Finally, we obtain the following three results on the existence of families of
periodic solutions for system (5).

######  Proposition 5.

There exists $K_{1}$ such that for any $k_{1}>K_{1}$ there exists
$\varepsilon=\varepsilon(k_{1})$ and for all $|F|<\varepsilon$ system (5) has
a two-parameter family of periodic solutions.

######  Proposition 6.

Let $k_{1}>0$, then there exists $\varepsilon=\varepsilon(k_{1})$, such that
for all $|m_{1}l_{1}-m_{3}l_{3}|<\varepsilon$ system (5) has a two-parameter
family of periodic solutions.

######  Proposition 7.

Let $k_{1}=0$ and $F<0$, then there exists $\varepsilon=\varepsilon(k_{1})$
such that for all $|m_{1}l_{1}-m_{3}l_{3}|<\varepsilon$ system (5) has a two-
parameter family of periodic solutions.

In all three propositions the proofs are based on the existence of families of
solutions which intersect the plane $\varphi_{1}=0$ transversely for $F=0$ and
$m_{1}l_{1}-m_{3}l_{3}=0$. From Proposition 4 we obtain that these families of
solutions exist for small $|F|$ and $|m_{1}l_{1}-m_{3}l_{3}|$.

From the existence of a two-parametric family of periodic solutions it follows
that locally (in a neighborhood of these solutions) the system is integrable:
there are two local first integrals, which can be defined in various ways.

(a)

(b)

Figure 3: Non-periodic (a) and periodic (b) solutions of the system for $F\neq
0$ and $m_{1}l_{1}\neq m_{3}l_{3}$. The parameters of the system are the same
in the both cases, initial conditions vary.

A natural question of whether or not it is possible to continue these two
first integrals globally to the whole three-dimensional space with coordinates
$\varphi_{1},v_{1},v_{2}$ is expected to have a negative answer. It follows
from the numerical results which show that in a general case for some initial
conditions there exist non-periodic trajectories lying on a two-dimensional
surface (Fig. 3). Therefore, there are less than two independent first
integrals.

The mechanism for the appearance of such trajectories is as follows. Let for
$F=0$ we have a family of periodic solutions, if $|F|$ is sufficiently small,
then all these solutions still intersect the plane $\varphi_{1}=0$ in two
points. Therefore, we obtain a family of periodic solutions. When $F$ becomes
sufficiently large, some solutions of this family do not intersect the plane
$\varphi_{1}=0$ and may become non-periodic. This breaks the integrability (or
super-integrability) of the system. In Fig. 4 a family of solutions starting
at the same point is shown. The only parameter of system that varies from
solution to solution is the magnitude of force $F$. For small $F$ all
solutions intersect the plane $\varphi_{1}=0$, but a non-periodic solution
appears when $F$ becomes larger than some critical value. This solution does
not intersect the considered plane.

Note that the above conditions for the existence of periodic solutions for
system (5) (for $F=0$) and (7) can be weakened and the families of periodic
solutions are much wider. We present simple conditions that can be obtained
without cumbersome calculations and detailed estimates.

Figure 4: The mechanism for the appearance of a non-periodic solution.

For the integrability in a four-dimensional space with an invariant measure,
the most natural definition of the integrability stems from the Jacobi
theorem: the system is integrable if it has two independent first integrals.

For any periodic solution in the reduced space with coordinates
$\varphi_{1},v_{1},v_{2}$, from the independent equation

$\dot{\varphi}_{2}=v_{2},$

it follows that in the four-dimensional space
$\varphi_{1},\varphi_{2},v_{1},v_{2}$ this solution lies on a two-dimensional
surface (the angular variables are considered modulo $2\pi$). If the solution
is not periodic in the reduced space $\varphi_{1},v_{1},v_{2}$, then its
projection from the four-dimensional space into a three-dimensional space
with, for instance, coordinates $\varphi_{1},\varphi_{2},v_{1}$ can densely
fill three-dimensional domains.

In the first case we have two first integrals. In the latter case, there may
be only one or no first integrals. As an illustration (Fig. 5) we present the
section of projections of two solutions into the coordinate space
$\varphi_{1},\varphi_{2},v_{1}$ by the plane $\varphi_{2}=\pi$: in the first
case the solution is periodic in the space with coordinates
$\varphi_{1},v_{1},v_{2}$, in the second case the solution is not periodic in
this space.

(a)

(b)

Figure 5: Regular (a) and irregular (b) behavior of solutions. The parameters
of the system are the same in the both cases, initial conditions vary.

In the conclusion of the section, we consider a sufficient generalization of
system (5) which is integrable. Indeed, the key property which is used in the
proves of integrability is the symmetry of the original system. Hence, if we
consider perturbations relatively small and satisfying the properties of
symmetry. For instance, let us consider the following system (here we assume
that $F=0$)

$\displaystyle\begin{split}&\dot{\varphi}_{1}=v_{1},\\\
&\dot{v}_{1}=\left(r_{1}-\frac{A_{12}}{A_{22}}r_{2}\right)\left(A_{11}-\frac{A_{12}A_{21}}{A_{22}}\right)^{-1}+f_{1}(\varphi_{1},v_{1},v_{2}),\\\
&\dot{v}_{2}=\left(r_{2}-\frac{A_{21}}{A_{11}}r_{1}\right)\left(A_{22}-\frac{A_{12}A_{21}}{A_{11}}\right)^{-1}+f_{2}(\varphi_{1},v_{1},v_{2}),\end{split}$
(8)

where $f_{i}(-\varphi_{1},v_{1},v_{2})=-f_{i}(\varphi_{1},v_{1},v_{2})$,
$i=1,2$ are sufficiently small. Then this system is integrable in the above
sense. For instance, a family of periodic solutions exists for system (8) if
we put $f_{1}=-\alpha\varphi_{1}\sin(v_{2})$,
$f_{2}=\alpha\varphi_{1}\sin(v_{1})$, for sufficiently small values of
$\alpha$. At the same time, there may exist asymptotically stable equilibria
in this system. Therefore, there is no smooth invariant measure in this system
and there cannot exist two first integrals, which are defined globally in the
whole phase space.

##  5 Case $k_{2}\neq 0$

In the above considerations we assume that the stiffness of the spring located
in the fixed point $O$ equals zero. This allows us to study the reduced system
with the three-dimensional phase space $\mathbb{S}\times\mathbb{R}^{2}$.

When $k_{2}\neq 0$ we have to consider the system defined on
$\mathbb{T}^{2}\times\mathbb{R}^{2}$, and this system has an invariant
measure. Based on the results from KAM theory for the reversible systems (for
details, please, see [11]), one can expect the co-existence of regions with
regular and chaotic behavior of the solutions.

A detailed numerical study of this case is out of the scope of this paper. We
only present some results showing that for the Ziegler pendulum the regular
dynamics can co-exists with the irregular behavior of trajectories.

After the change of time variable

$t\mapsto t\cdot(A_{11}A_{22}-A_{12}A_{21})$

system (5) can be presented in the following form (here we also the equation
for $\dot{\varphi}_{2}$)

$\displaystyle\begin{split}&\dot{\varphi}_{1}=v_{1}(A_{11}A_{22}-A_{12}A_{21}),\\\
&\dot{\varphi}_{2}=v_{2}(A_{11}A_{22}-A_{12}A_{21}),\\\
&\dot{v}_{1}=A_{22}r_{1}-A_{12}r_{2},\\\
&\dot{v}_{2}=A_{11}r_{2}-A_{12}r_{1},\end{split}$ (9)

where we use the same notations as above.

One of the most popular methods of numerical study of chaos is $mLCE$ (a
simple method based on the estimation of the maximum Lyapunov characteristic
exponent). A detailed explanation of this approach can be found in [12, 13,
14].

For a given solution, the Lyapunov characteristic exponents can be calculated
as follows

$\chi(X(t),x)=\limsup\limits_{t\to+\infty}\frac{1}{t}\ln\frac{\|X(t)x(0)\|}{\|x(0)\|},$

where $X(t)$ is the fundamental matrix of the linearized system for system
(9). The system is linearized along the considered solution. The evolution of
vector $x(t)$ is described by the following linear system

$\dot{x}=A(t)x.$

For our system, we have

$A(t)=\begin{pmatrix}v_{1}A_{11}\displaystyle{\frac{\partial
A_{22}}{\partial\varphi_{1}}}-2v_{1}A_{12}\frac{\partial
A_{12}}{\partial\varphi_{1}}&0&A_{11}A_{22}-A_{12}A_{21}&0\\\
v_{2}A_{11}\displaystyle{\frac{\partial
A_{22}}{\partial\varphi_{1}}}-2v_{2}A_{12}\frac{\partial
A_{12}}{\partial\varphi_{1}}&0&0&A_{11}A_{22}-A_{12}A_{21}\\\
\displaystyle{\frac{\partial
A_{22}}{\partial\varphi_{1}}r_{1}+A_{22}\frac{\partial
r_{1}}{\partial\varphi_{1}}-\frac{\partial
A_{12}}{\partial\varphi_{1}}r_{2}-A_{12}\frac{\partial
r_{2}}{\partial\varphi_{1}}}&-A_{12}\displaystyle{\frac{\partial
r_{2}}{\partial\varphi_{2}}}&A_{22}\displaystyle{\frac{\partial
r_{1}}{\partial v_{1}}}-A_{12}\displaystyle{\frac{\partial r_{2}}{\partial
v_{1}}}&A_{22}\displaystyle{\frac{\partial r_{1}}{\partial
v_{1}}}-A_{12}\displaystyle{\frac{\partial r_{2}}{\partial v_{1}}}\\\
\displaystyle{\frac{\partial
A_{22}}{\partial\varphi_{1}}r_{1}+A_{22}\frac{\partial
r_{1}}{\partial\varphi_{1}}-\frac{\partial
A_{12}}{\partial\varphi_{1}}r_{2}-A_{12}\frac{\partial
r_{2}}{\partial\varphi_{1}}}&A_{11}\displaystyle{\frac{\partial
r_{2}}{\partial\varphi_{2}}}&A_{22}\displaystyle{\frac{\partial
r_{1}}{\partial v_{1}}}-A_{12}\displaystyle{\frac{\partial r_{2}}{\partial
v_{1}}}&A_{22}\displaystyle{\frac{\partial r_{1}}{\partial
v_{1}}}-A_{12}\displaystyle{\frac{\partial r_{2}}{\partial
v_{1}}}\end{pmatrix}$

In practice, the following equation is used to obtain the maximum Lyapunov
characteristic exponent

$\chi(X(t))=\lim\limits_{t\to+\infty}\frac{1}{t}\ln\frac{\|X(t)x(0)\|}{\|x(0)\|}.$

It is known that for almost all initial vectors $x(0)$ and for a wide range of
matrices $X(t)$ this equation gives the maximum Lyapunov characteristic
exponent.

Let us consider the system with the following values of mechanical parameters:

$\begin{split}&m_{1}=1,\\\ &m_{2}=1,\\\
&m_{3}=3/2,\end{split}\quad\begin{split}&l_{1}=1,\\\ &l_{2}=1,\\\
&l_{3}=1,\end{split}\quad\begin{split}&k_{1}=1,\\\ &k_{2}=1,\\\
&F=2,\end{split}$

and two solutions with initial conditions $(\pi,0,1/10,1/10)$ and
$(\pi,0,20,20)$ (the order of the coordinates is as follows:
$\varphi_{1},\varphi_{2},v_{1},v_{2}$). Below we present the results of
calculation of $mLCE$ for these solutions (Fig. 5). It can be seen that in one
case the value of $mLCE$ always remains positive (a property of a chaotic
trajectory) and for another solution this value tends to zero.

(a) $mLCE=0$

(b) $mLCE>0$

Figure 6: Zero and positive maximum Lyapunov characteristic exponents.

Regularity and irregularity of solutions can be also seen in the projection of
these solutions into the three-dimensional space $(\varphi_{1},v_{1},v_{2})$.
In the case of regular behavior with zero $mLCE$, the solution lies on a two-
dimensional surface. For the trajectory with strictly positive $mLCE$, th
closure of the projection of the trajectory covers three-dimensional domains
(Fig. 6).

(a)

(b)

Figure 7: Regular (a) and irregular (b) behavior of the solutions. The
parameters of the system are the same in the both cases, initial conditions
vary.

We also present ‘the Poincare sections’ of a family of solutions, from which
it is also can be seen that the regular dynamics co-exists with the chaotic
dynamics. We consider ten initial conditions $(\pi,0,1/10,(n+1)/10)$,
$0\leqslant n<10$, $n\in\mathbb{Z}$. For each corresponding solution we
consider its projection into the three-dimensional space
$(\varphi_{1},v_{1},v_{2})$, and then the section of this projection by the
plane $v_{1}=0$.

Figure 8: Regular and chaotic solutions.

The obtained sets of points are not the Poincare section in the usual sense
since we consider the section of the projections of solutions. In particular,
it is possible that there exist more than one solution which intersects the
plane $v_{1}=0$ in a given point. The resulting sections are presented in Fig.
7. It can be seen that for relatively small initial velocity the solutions are
regular (lie on a two-dimensional surface) and become chaotic as velocity
becomes larger.

##  6 Conclusion

The dynamics of a planar double pendulum in the presence of a follower force
is considered. The system is studied in the case when the governing system of
equation can be reduced to a system in a three-dimensional space, and in the
general case when the system is considered in a four-dimensional space. The
following results are obtained:

  1. 1.

It is proved that, in the case of zero stiffness of the spring located in the
point of suspension, there exist two-parametric families of periodic solutions
(for some values of the mechanical parameters).

  2. 2.

It is shown numerically that the periodic solutions do not fill the whole
phase space. In particular, some solutions lie on two-dimensional surfaces and
the system does not have two independent first integrals defined globally in
the space with coordinates $\varphi_{1},v_{1},v_{2}$. Moreover, the system is
not integrable in the sense of the Jacobi theorem if we consider the dynamics
in the original (not reduced) phase space
$\mathbb{T}^{2}\times\mathbb{R}^{2}$.

  3. 3.

It is shown that in the general case of the Ziegler pendulum when $k_{2}\neq
0$ the regular and the chaotic dynamics co-exist in different regions of the
phase space.

Let us briefly discuss the above results and outline some possible directions
of further research. First, let us note that we make the conclusion on the
integrability of the system without explicitly presenting the required set of
the first integrals. We only prove that there exist periodic solutions for the
system and, therefore, the system can be locally integrable. Moreover, there
is a strong evidence that these local first integrals cannot be continued to
the whole phase space. It might be interesting to understand this result
deeper: to the best of our knowledge there are no other mechanical systems
that are integrable in one part of the phase space and not integrable, or
integrable in another sense, in the rest of the phase space.

An important property of the circulatory system (5), which cannot be seen in
non-integrable Hamiltonian systems, is that there are regions in the phase
space of system (5) such that all invariant tori are not destroyed when we add
a perturbation to the system. For Hamiltonian systems it is typical that
resonant tori are destroyed by perturbations. Since resonant tori are dense in
the phase space, we can usually conclude that the system becomes non-
integrable.

At the same time, the co-existence of the regular and the chaotic behavior of
solutions is typical for Hamiltonian systems (see, for instance, [15, 16]).
Nevertheless, we have shown that chaotic trajectories in circulatory systems
(the Ziegler pendulum for $k_{2}=0$) may appear in a different manner
comparing to Hamiltonian systems. The existence of chaotic regions in analytic
Hamiltonian systems means that the system is not integrable, yet the existence
of local set of independent first integrals is still possible if we have
irregular behavior of solutions in a circulatory system.

## Acknowledgments

The author is grateful to Valery Kozlov for his invaluable advice and useful
suggestions.

The work was supported by the Russian Science Foundation (Project No.
19-71-30012).

## References

  * [1] Hans Ziegler  “Die stabilitätskriterien der elastomechanik”  In _Ingenieur-Archiv_ 20.1 Springer-Verlag, 1952, pp. 49–56
  * [2] G Herrmann and RW Bungay  “On the stability of elastic systems subjected to nonconservative forces”, 1964
  * [3] VV Bolotin and NI Zhinzher  “Effects of damping on stability of elastic systems subjected to nonconservative forces”  In _International Journal of Solids and Structures_ 5.9 Elsevier, 1969, pp. 965–989
  * [4] Oleg N Kirillov  “Nonconservative stability problems of modern physics”  In _Nonconservative Stability Problems of Modern Physics_ de Gruyter, 2013
  * [5] Troy Shinbrot, Celso Grebogi, Jack Wisdom and James A Yorke  “Chaos in a double pendulum”  In _American Journal of Physics_ 60.6 American Association of Physics Teachers, 1992, pp. 491–499
  * [6] Tomasz Stachowiak and Toshio Okada  “A numerical analysis of chaos in the double pendulum”  In _Chaos, Solitons & Fractals_ 29.2 Elsevier, 2006, pp. 417–422
  * [7] Holger R Dullin  “Melnikov’s method applied to the double pendulum”  In _Zeitschrift für Physik B Condensed Matter_ 93.4 Springer, 1994, pp. 521–528
  * [8] Valery V Kozlov  “On the Integrability of Circulatory Systems”  In _Regular and Chaotic Dynamics_ 27.1 Springer, 2022, pp. 11–17
  * [9] Valery V Kozlov  “The Euler-Jacobi-Lie integrability theorem”  In _Regular and Chaotic Dynamics_ 18.4 Springer, 2013, pp. 329–343
  * [10] Valery V Kozlov  “Integrals of Circulatory Systems Which are Quadratic in Momenta”  In _Regular and Chaotic Dynamics_ 26.6 Springer, 2021, pp. 647–657
  * [11] Mikhail B Sevryuk  “Reversible systems”  Springer, 2006
  * [12] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli and Jean-Marie Strelcyn  “Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory”  In _Meccanica_ 15.1 Springer, 1980, pp. 9–20
  * [13] Allan J Lichtenberg and Michael A Lieberman  “Regular and chaotic dynamics”  Springer Science & Business Media, 2013
  * [14] Ch Skokos  “The Lyapunov characteristic exponents and their computation”  In _Dynamics of Small Solar System Bodies and Exoplanets_ Springer, 2010, pp. 63–135
  * [15] Vladimir I Arnold, Valery V Kozlov and Anatoly I Neishtadt  “Mathematical aspects of classical and celestial mechanics”  Springer Science & Business Media, 2007
  * [16] Valerij V Kozlov  “Symmetries, topology and resonances in Hamiltonian mechanics”  Springer Science & Business Media, 2012