CERN antiproton target: hydrocode analysis of its core material dynamic response under proton beam impact

Antiprotons are produced at CERN by colliding a 26 GeV/c proton beam with a fixed target made of a 3 mm diameter, 55 mm length iridium core. The inherent characteristics of antiproton production involve extremely high energy depositions inside the target when impacted by each primary proton beam, making it one of the most dynamically demanding among high energy solid targets in the world, with a rise temperature above 2000 {\deg}C after each pulse impact and successive dynamic pressure waves of the order of GPa's. An optimized redesign of the current target is foreseen for the next 20 years of operation. As a first step in the design procedure, this numerical study delves into the fundamental phenomena present in the target material core under proton pulse impact and subsequent pressure wave propagation by the use of hydrocodes. Three major phenomena have been identified, (i) the dominance of a high frequency radial wave which produces destructive compressive-to-tensile pressure response (ii) The existence of end-of-pulse tensile waves and its relevance on the overall response (iii) A reduction of 40% in tensile pressure could be obtained by the use of a high density tantalum cladding.


I. INTRODUCTION
The AD-Target system is the main particle production element of the CERN Antiproton Decelerator. Antiprotons are produced by colliding a proton beam of 26 GeV/c from CERN Proton Synchrotron (PS) with a fixed target made of a dense and high-Z material. The high energy collision of the proton beam with the atoms of the target creates a shower of secondary particles, and among them, antiprotons. The generated antiprotons travel through the target assembly and are collected downstream by a magnetic focusing device using their charge, after which, they are selected and magnetically conducted via the injection line to the Antiproton Decelerator (AD) Complex for antimatter research experiments. Antiparticles should be identical to matter particles except for the sign of their electric charge. It is well known, however, that matter and antimatter are not exact opposites; nature seems to have a one-part in 10 billion preference for matter over antimatter. Understanding matterantimatter asymmetry is one of the greatest challenges in physics today.
The relevance of the antiproton target is not only motivated by the significance of its purpose -to produce new particles to serve as a window to new physics-but for the challenge that involves to overcome the engineering problem of its design. The characteristics of antiproton production require a very compact target in order to avoid antiproton reabsorption in the surrounding material and to be as close as possible to a punctual source for the antiproton downstream collector system [1,2]. For this reason, a very thin rod of a high density material and a very focused primary proton beam have to be used. This results in extremely high energy depositions reached inside the target core as a consequence of each proton beam impact, which makes the antiproton target a unique case even among high energy particle-producing primary targets. As a matter of fact, calculations show that in each proton pulse impact the peak temperature rises up to 2000 • C. Dynamic pressure waves of several GPa and strain rates well above 10 4 s −1 are easily exceeded as a consequence of the sudden expansion of the material. These waves propagate through the target and surrounding containment, where expected stresses are also considerably high [3]. All these phenomena turn into a real challenge to predict the target material response and behavior during operation. In addition, radiation damage including material embrittlement, swelling and noble gas production inside the target material may also play an additional role [4][5][6].
These damaging processes, however, are not considered in this publication and are subjects The present configuration of the antiproton target dates back to the late 80's and it was obtained after more than 10 years of iterations and improvements of the material choice and conceptual designs. Copper, tungsten, rhenium and iridium were employed during this decade. The current design consists of a core target in form of a 55 mm long, 3 mm diameter iridium rod -one of the densest elements on earth, second only to osmium-encapsulated in a 15 mm diameter graphite matrix; the whole assembly is in turn embedded into a watercooled Ti-6Al-4V body as shown in Figure 1. After 2016 a new upgrade of the AD facility is planned, including an optimized re-design of the full antiproton production system in order to enhance its performance, reliability as well as to guarantee the continuity of the antiproton physics for the next 20 years [7][8][9].
The study of the dynamic response of the target material is fundamental in order to reduce the uncertainties concerning the target core structural state and to provide a robust design. Any damage and effective loss of density of the target material will lead to target exchange since the proton-target interactions and antiproton production will be reduced to the point that the target has to be replaced periodically to maintain the desired performance [10,11]. Unfortunately, post examination of failed targets is very complex due to high activation of the irradiated material. For this reason, a reliable methodology based on advanced computational tools and experimental tests is proposed in order to reduce uncertainties on the target response and to assess the choice of future design and candidate material. In the present study, the numerical part of this approach is shown with the aim of identifying the AD-target material working conditions and analyze its governing phenomena.
The thermo-mechanical effects induced in materials by high energy particle beam impact 3 have been extensively studied in accelerator technology in a wide number of facilities around the world for the design of BIDs (beam intercepting devices) such as targets, collimators and dumps. References [12][13][14] are a few among many examples. In these studies, the energy deposition map arising as consequence of proton beam interaction with matter is calculated by means of Monte Carlo simulations and then applied as an internal energy load to thermo-mechanical FE solvers. The thermo-mechanical effects induced in solid were also studied analytically in [15], where dynamic stresses were solved for a single circular traversal section assuming impact of a Gaussian beam. All these calculations, however, are normally performed under assumptions of material elastic regime and considering quasi-static analysis, focused towards steady state conditions since raise of temperatures are within the order of hundreds of degrees per pulse. Unfortunately, these assumptions are by far out of reality for our problem of study, where deposited energy is much higher and fast. The material will certainly deform plastically or break. As a matter of fact, only the target which was operating for antiproton production at US Fermilab has worked in a similar range of energy deposition density as CERNs AD-target [16].
Hydrocodes are a family of advanced highly non-linear computational tools which are able to fully simulate the dynamic response of materials [17]. These codes numerically solve the equations arising from mass, momentum and energy conservation. They are able to simulate the material response beyond plasticity and fracture taking into account strain rate and temperature dependence by coupling these equations with material constitutive laws as equations of state (EOS), strength and failure models. Hydrocodes have been historically used in ballistic and military applications where the dynamic load originates from high velocity impacts and explosions. However, they have recently started being applied in accelerator technology for simulations of BIDs, where the load is generated by the sudden thermal expansion of the material impacted by particle pulses and the subsequent propagation of the produced stress waves. Current programs for development of accelerator technologies require higher and higher particle beam energies and intensities, which push engineers to use these codes to investigate mechanical damage in materials hit by high intensity beams as a consequence of abnormal operation, in which hydrodynamic behaviour of materials is expected. A comprehensive analysis using AUTODYN R was done for the study of the structural behaviour of the Main LHC absorber block in case of a total beam dilution failure [18]. Numerical studies of high dynamic transient effects of pulse beams on BIDs can be found in [19,20], where hydrodynamical calculations for uranium beam impacted on copper targets for ISOLDE were performed. Very recent publications for research on new collimator materials using ANSYS AUTODYN R can be found in [21,22], while using as well LS-DYNA R in references [23][24][25]. In the same way, literature including hydrodynamic calculations using BIG2 code for the design of solids targets for the FAIR facility in Germany can be found in [26,27].
The goal of the present paper is therefore, the application of hydrocodes in order to simulate the dynamic response of the AD-Target core after one proton pulse impact, in order to gain important information concerning its response, understand the dominating phenomena and learn lessons which could define strategies to provide a new design which should be as robust as possible. Figure 2 shows a scheme of the methodology employed in this study. First, an energy deposition map in the target material as a consequence of the proton beam-target atoms interaction is calculated by means of FLUKA Monte Carlo simulations, a highly advanced particle transport code developed at CERN which is able to simulate the energy released by these nuclear reactions also taking into account secondary reactions of the particle cascade generated by the proton beam when going through the target material [28]. This energy deposition map is then applied as an internal energy load to the commercial hydrocode ANSYS AUTODYN R via an internal user subroutine and the user function EXSIE3. The energy is applied consistently with the PS beam parameters, whose proton beam is composed by 4 proton bunches spaced by 105 ns, bunch length of 30 ns for a total pulse intensity of 1.5·10 13 protons per pulse [9]. The total energy deposited in the target core is approximately 1.34 kJ, which results in 11.17 GW considering that it is deposited only in 120 ns. Moreover, the power is applied in a small volume of the target core, producing a mean power density of 2.87 · 10 10 T W/m 3 .

II. COMPUTATIONAL APPROACH
At this point, it is important to note that the interaction between FLUKA and ANSYS AUTODYN is only a one-way coupling, which means that changes in the deposited energy during the pulse impact as a consequence of the change of target material density are not taken into account.   [17]. Unfortunately, the material constants associated to many of these models are rarely available in literature, and are often classified as it is drawn from military research. Moreover, most of the AD-target material candidates such as iridium or tungsten-rhenium alloys are uncommon and they have never been studied at any dynamic regime at all, even less so for the ones present in the case of study, with temperatures above 2000 • C and strain rates above 10 4 s −1 . It is for this reason that the calculations here are done assuming pure tungsten as a target material, for which more material models are available.
Nevertheless, tungsten was and is a strong candidate material and the lessons learned from this study depend mainly on the pure dynamic response, making it possible to extrapolate theses conclusions to other candidate materials. Furthermore, experimental tests under real proton beam and similar conditions as reached in the AD-target are also foreseen to validate these phenomena.

Equation of State (EOS)
An equation of state is a constitutive law which expresses the relation between thermodynamic state variables as pressure, energy and density. EOS were originally developed for ideal gases and subsequently extended to all states of matter [29]. In continuum mechanics, the stress tensor is usually decomposed in a hydrostatic component and deviatoric component. The EOS governs hydrostatic component associated with the change of pressure as a function of volume and internal energy and it is strongly correlated with the wave propagation velocity and shock response. Ideally, the EOS includes also change of phase of the material, although in practice it is difficult to include them in a single formulation. There are different EOS formulations depending on the regime of applications and assumptions, Linear EOS, Mie-Grüneisen, Puff, Tillotson EOS [30] and SESAME tables [31] are some of them. For this study Mie-Grüneisen is employed. This model is one of the most widely used EOS and one of which more parameters are available in the literature, working especially well in the regime of not very high pressures where the relation between shock speed and particle speed can still be considered linear. The AUTODYN R implementation of Mie-Grüneisen takes as well into account differences between compression and tensile states.
This formulation does not take into account change of phase, which is perfectly valid for the case under study -where melting of material does not occur due to the high melting point of the candidate materials-and the pressures reached are in the order of GPa's, which in terms of shock physics, are relatively low.

Strength Model
The strength model, on the other hand, defines the deviatoric part of the stress-strain tensor, i.e. it governs the changes in shape as well as deformation beyond plasticity (flow stress). The strength model is therefore a constitutive relation that links stress with strain for different strain rates and temperatures. There are a wide number of strength model formulations depending on strain rate-temperature influence assumptions, type of material, and regime of application. Some of them are purely empirical as Johnson-Cook [32], which is one of the simplest and more widely used strength model, or more physical based as Zerillini-Amstrong [33] which is based on dislocation dynamics. Other formulations are Steinberg-Guinan [34], Steinberg-Lund [35] or the mechanical threshold model [36]. For the present study, Johnson-Cook strength model was used, which takes into account the temperature and strain rate dependence on the material as well as its response beyond plasticity. This model is relatively simple, since it considers the influence of temperature and strain rate in an uncoupled way, but has been shown to provide good agreement between experiments and simulations in other studies of proton beam impact phenomena on materials [21]. The parameters employed for the J-C strength model of pure tungsten were obtained by Hopkinson bar tests at different high temperatures and high strain rates by K.T Ramesh et al. [37] and are partially extrapolated by the code to the regime reached in the AD-target.

Failure Model
Finally, the third necessary element to simulate realistic material behavior is a failure model. A material cannot withstand local tensile stresses greater than its limits. When this threshold is exceeded, the code assumes that the material can no longer sustain any shear stress or any tensile pressures in the corresponding element. There are two main types of failure models used by the hydrocodes; (i) the material will instantaneously fail when it locally overcomes a limit for one or a certain number of variables (such as strain to fracture, tensile hydrostatic stress, maximum principal stress) and (ii) The failure model is based on the cumulative damage of the material when a certain variable is re-iteratively exceeded [38]. The first one is generally associated to brittle failure mode while the second with ductile fracture [43]. In the present study, a model corresponding to the first type is employed, named the Minimum Hydrostatic Pressure model [39]. This model assumes that Simulations are performed using the reference geometry shown in Figure 3 (a), which corresponds to the target core and containing graphite matrix. In addition to this configuration, two other configurations have been investigated, including a copper and tantalum cladding up to 1 mm thickness around the target core, (Figure 3    than 800 MPa, well exceeded in this study. Therefore, the pressure values shown in these plots can only be taken into account as a qualitative way since the energy diffused by plastic deformation processes is not taken into account. Simulations assuming elastic material, however, allow a better understanding of the pressure wave propagation as no damping of the waves on the material occurs. Figure 5 shows the pressure response in the target during the first 7 µs . A four stepped rise of pressure corresponding to the four impacts of four proton bunches spaced by 120 ns can be observed. Then, a clear high frequency pressure wave response takes place. This wave has a period of 0.85 µs , corresponding to the time needed for the induced expansion to travel and come back from the rod circumferential surface. The pressure wave in the centre and periphery (gauges 7 and 9) are practically in phase, meaning that the whole section of the rod expands simultaneously. This is due to the fact that the sudden rise of temperature in the material does not only happen at the center but also at the periphery of the rod, where the increase of temperature is also significant (above 300 • C, as can be seen in figure 4). The consequence is that the pressure waves are not just traveling from a central, localized area, but there are infinitesimal pressure waves with infinitesimal departing points along the rod material. The amplitudes of these pressure waves are attenuated as they travel to the peripheral surface, and amplified as they do it through the center due to the cylindrical geometry. The pressure, and other variables, oscillations are the result of the interactions of all these waves at each single point. rod is therefore longitudinally expanding in a simultaneous way within all its volume, in a similar manner as it was happening radially. The reason for this phenomenon, again, is that the temperature rise caused by the deposited energy is also covering a great part of the longitudinal volume of the core rod. An important point to be detailed here is that the velocity of propagation of the radial and the longitudinal waves do not match. This can be easily found out by just dividing the longitudinal and radial period by its corresponding radial and longitudinal distances. Initially one may think that the nature of the origin of the wave is the tri-axial expansion of the inner material and therefore the wave should travel through the bulk material at the same speed, no matter the direction. This however is not perfectly true, since ultimately, the tri-axial expansion in a single point and its propagation is governed and constrained by the surrounding material as well. In that sense, and in this specific case where the geometry is that small, the amount of surrounding material does matter. In fact, the resulting radial oscillation has a faster speed of propagation than the longitudinal one since there is much less surrounding material limiting the expansion radially meaning that the rod material simultaneously expands and contracts within its transversal section than longitudinally. The previous subsection allowed pointing out the existence of an end-of-bunch tensile wave taking place due to the inertia effects. We have investigated whether these tensile waves can have an important influence on the overall response. For this reason, a parametric analysis assuming proton impacts with different pulse length has been carried out. For this analysis it is assumed that all the energy is deposited in the target in continuous pulses (single bunch). The total energy applied for each of the scenarios is the same, with the only difference of the duration of the pulse. A Johnson-Cook strength model including plastic deformation of the material is used. Figure 10 shows the pressure response in the center of the target assuming times of proton pulses from 0.1 µs (blue) to 1.2 µs (pink) in steps of 0.2 µs. As can be seen in the plot, the magnitude of pressures reached change completely depending on the length of the pulse, ranging from +/-9 GPa with 0.1 µs pulse to +1 GPa 16 only in compressive stress with 1 µs . The reason for this change is not only the evident fact that the shorter the pulse the greater the power and inertia but the influence of when the pulse is finishing. It can be observed in the plot that a tensile wave which changes the trend of the pressure wave always takes place at the end of the pulse. The effect of this end-ofpulse tensile wave differs significantly depending on whether the material is going to tensile or compressive states. When the material is going to tensile states, the end-of-pulse wave and radial wave produce a constructive interference leading to very large tensile stresses in the first period of the radial wave, which propagates to the rest of the transient pressure response. This can be easily appreciated in figure 11; higher end-of-pulses pressure drops take place in the cases where pulses finish when DP/Dt < 0. This is the case of pulses of 0.43 µs length (the real of the AD-target), 0.6 µs length and 1.2 µs length. On the other hand, in the pulses which end between 1/2 and 3/4 of the period of the radial wave (when the pressure is going to compressive states and therefore DP/Dt > 0) the end-of-pulse tensile wave is counteracted by the radial wave. This is the case of 0.8 µs length and 1 µs pulses in figures 10 and 11, where the tensile pressure reached is below -1 GPa. This phenomenon is even clearer when comparing the 1 µs length pulse (yellow) with the one of 1.2 µs. Higher tensile and compressive stresses are reached in the latter, even if the energy is deposited slower in this case. The only exception to this response takes place with the pulse 0.1 µs, where due to the high inertia, the amplification of the end-of-pulse tensile wave takes place even if DP/Dt > 0. In any case, a straightforward conclusion that can be drawn from this analysis is that a way to reduce the harmful tensile stresses reached in the AD-target, given that the length of the pulses is strongly constrained by the PS and AD operation, would be to increase the diameter of the target core to 5-6 mm in order to increase the period of the radial wave and avoid its constructive interference with the end-of-pulse target. This diameter increase is however limited by the antiproton re-adsorption at the periphery of the target, fact which demands detailed trade-off studies on this possibility.

D. Influence of Cladding and application of Failure Models to the Cases of Study
Another studied strategy in order to decrease the magnitude of the tensile pressure in the center of the target core could be the addition of an external cladding. It is well known in the shock and explosive physics literature that the use of multi-layer armatures can effectively 0 0. 2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6  reduce the shock-induced tensile stresses in explosively expanding materials. The key point resides in the acoustic impedance mismatch between the different materials in a multi-layered media, which governs the fraction of the pressure wave energy that is transited/reflected between them [42,43]. The greater the impedance mismatch, the greater the portion of the shock energy reflected at the boundary, and the smaller, the greater portion transmitted to the adjacent media. The ideal situation from the shock wave view point would be that no impedance mismatching exist between the target core and surrounding material, therefore all the energy of the radial pressure wave is transmitted to the adjacent media instead of being reflected to the core as a tensile and destructive wave. This possibility is studied numerically in the present section by assuming three different scenarios: • No cladding between the 3 mm diameter high density target core and the surrounding graphite matrix (Figure 3-a). This is the case of the current design of the target, where the impedance mismatch is very big due to the large difference between tungsten and graphite's density and sound velocity (≈ 30 times ratio in acoustic impedance between target core and graphite) • Copper cladding of 1 mm thickness (Figure 3  Assuming Jonson-Cook strength model for the tungsten core and models considering plasticity in the cladding consistently with table I. No failure models are considered. The plot clearly shows a reduction of tensile pressure of 44% when using Ta cladding pressure in the radial wave is significantly reduced when using copper and tantalum cladding, by 44 % when using tantalum. The reason for this reduction is probably a combination of 3 phenomena: first, the fact of the better impedance match between the core of tungsten and tantalum which causes a higher fraction of the energy of the pressure wave to be transmitted to the cladding instead of coming back as a tensile wave. Second, the shift in the period of the radial wave due to the larger diameter as a consequence of the cladding, which could partially avoid the constructive interference between the end-of-pulse wave and the radial ones explained in the previous section. Third, the plastic deformation in the tantalum or copper cladding material, which absorbs energy from the pressure wave, avoiding its reflexion to the target core. This reduction by 44 % in the maximum tensile pressure reached in the center of the target when using tantalum can be very important from the point of view of target survival and increase of antiproton yield, especially if it is enough to reduce the tensile pressure below the spall strength of the material (which in tungsten was identified at 2.6 GPa). This can be easily seen in figure 13, where minimum hydrostatic pressure failure model is applied in the three cases of interest after a single pulse of proton beam. As observed in the figure, for the case of a tungsten target without cladding, a large part of the volume of the target core is fragmented, leading to a reduction of antiproton production due to loss of effective 20

IV. CONCLUSIONS
In the present study a detailed analysis of CERN antriproton decelerator target has been carried out. This analysis includes a summary of the existing literature of the current design, which dates from the late 80s, as well as of the previous designs built during the 10 years of iterations which led to it. This study emphasizes the extreme conditions reached in the 3 mm diameter, high density target core, in terms of adiabatic increase of temperature of more than 2000 • C each time the proton beam hits the target and the subsequent pressure wave and dynamic response. These facts probably make the AD-target the most demanding high energy target currently in operation around the world.

21
The understanding of the dynamic phenomena induced on the target core as a consequence of proton beam impact turns out to be fundamental in order to provide a robust future design for the next 20-30 years of operation of the antiproton decelerator. Hydrocodes are a very powerful numerical tool, historically used for shock and impact physics, which are the only way to properly simulate this dynamic response of the target. The application of them to the target core is even more interesting, since its small and simple cylinder geometry allows to easily appreciate the propagation of the generated pressure waves and understanding fundamental aspects about them. Simulations starting from simple models assuming perfectly elastic material, while adding gradually more complex models which take into account plasticity or fracture, helped as well for this purpose. In this way, three main phenomena have been identified and described in this study: (i) The identification of radial and longitudinal pressure waves present in the target core with a period of 0.85 µs and 24 µs respectively, the high frequency radial wave being the most important and critical one since the longitudinal one is rapidly damped by the plastic deformation of the material. The radial wave on the other hand produces huge oscillating compressive-to-tensile pressure response in the center and periphery of the target core, probably leading to its internal fracture in the first oscillations, when reaching tensile stresses above the spall strength of the material.
(ii) The existence and importance of the end-of-pulse tensile waves, consequence of the inertia of the heated material expansion. These tensile waves become especially important when they are in phase with the tensile part of the compressive-to-tensile radial wave described above. This is the case of the AD-target core due to its small diameter of only 3 mm.
Avoiding this constructive interference will be a way to reduce this limiting phenomenon.
(iii) The possibility of using a 1 mm thickness cladding surrounding the 3 mm target core as a way to decrease the acoustic impedance mismatch between the target core and graphite containing matrix, with the aim of decreasing the tensile pressure reached in the center of the high density target core. Simulations show that a reduction of up to 40 % in the maximum tensile pressure can be achieved by a tantalum cladding.
Nevertheless, one must be aware of the intrinsic limitations of these simulations, which are using material models based on simple assumptions and experimentally fitted as Jonhson-Cook can be. These models are applied here to a much more complex case than a simple mechanical impact test. The nature of the load itself, the variation of temperature and strain rates in this short time, will lead to lots of phenomena which cannot be considered by the code, like possible re-crystallization processes if there is enough time for it, presence of pre-existing defects, material mechanical hysteresis, radiation damage . . . as well as the fact that material strength models used are partially extrapolated since there is no availability in the literature of them at this extreme regime.
For this reason, the analysis in the current study only represents the numerical approach of the antiproton target re-design process. An experiment in the testing facility HiRadMat at CERN is foreseen to experimentally complete this numerical work. In this experiment different target materials as Ir, W, Mo, Ta, TZM will be irradiated with high energy proton beam from the SPS, recreating the same conditions as reached in the real AD-target and exposed here. The goal of this experiment will be to validate these calculations, gain experimental insights of the real material response in these conditions and assess new potential candidates target materials for a future design.

V. ACKNOWLEDGMENTS
The authors want to thank Herta Richter for the development of the subroutine used to import the energy deposition from FLUKA to AUTODYN R . In addition, the authors gratefully acknowledge the advices of Professor Lorenzo Peroni, which helped with the interpretation of some of the results, as well as the English proofreading of this document kindly done by Ayse Karatepe.