# TRIM Simulations Tool for $\mu^{+}$ Stopping Fraction in Hydrostatic
Pressure Cells

Frank Elson1,∗ Debarchan Das2 Gediminas Simutis3,4 Ola Kenji Forslund4 Ugne
Miniotaite1 Rasmus Palm1 Yasmine Sassa4 Jonas Weissenrieder1 and Martin
Månsson1,+ 1Department of Applied Physics, KTH Royal Institute of Technology,
SE-106 91 Stockholm, Sweden  2Laboratory for Muon Spin Spectroscopy, Paul
Scherrer Institut, CH-5232 Villigen PSI, Switzerland  3Laboratory for Neutron
and Muon Instrumentation, Paul Scherrer Institut, CH-5232 Villigen PSI,
Switzerland  4Department of Physics, Chalmers University of Technology,
Göteborg, SE-412 96, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

For quantum systems or materials, a common procedure for probing their
behaviour is to tune electronic/magnetic properties using external parameters,
e.g. temperature, magnetic field or pressure. Pressure application as an
external stimuli is a widely used tool, where the sample in question is
inserted into a pressure cell providing a hydrostatic pressure condition. Such
device causes some practical problems when using in Muon Spin
Rotation/Relaxation ($\mu^{+}$SR) experiments as a large proportion of the
muons will be implanted in the pressure cell rather than in the sample,
resulting in a higher background signal. This issue gets further amplified
when the temperature dependent response from the sample is much smaller than
that of the pressure cell,which may cause the sample response to be lost in
the background and cause difficulties in aligning the sample within the beam.
To tackle this issue, we have used pySRIM [1] to construct a practical and
helpful simulation tool for calculating muon stopping fractions, specifically
for the pressure cell setup at the $\mu$E1 beamline using the GPD spectrometer
at the Paul Scherrer Institute, with the use of TRIM simulations. The program
is used to estimate the number of muon stopping in both the sample and the
pressure cell at a given momentum. The simulation tool is programmed into a
GUI, making it accessible to user to approximate prior to their experiments at
GPD what fractions will belong to the sample and the pressure cell in their
fitting procedure.

##  1 Introduction

High pressure studies using muon spin resonance ($\mu^{+}$SR) has become very
popular as a ’clean’ approach for probing and tuning material properties, and
especially within the field of magnetic materials and quantum criticality [2,
3, 4]. The purpose of this is to investigate the phase diagram of a given
material as a function of an external perturbation (i.e. pressure but also
e.g. temperature and magnetic field). Through such an approach, it is possible
to shed light onto a plethora of complex phenomena, such as unconventional
superconductivity [5], quantum spin liquids [6] or frustrated magnetism [7].
The combination of pressure and the $\mu^{+}$SR technique, in which a
microscopic understanding of the ground state of these quantum materials can
be obtained, provides a uniquely powerful tool to get insight into the
intrinsic properties of the material. However, to apply hydrostatic pressure
(e.g. at GPD, which will be the spectrometer used as the example throughout
this paper), the samle needs to be compressed into a pellet and inserted into
the pressure cell made of either MP45N or Cu-Be. Before starting any
measurements, the muon beam and its momentum need to be adjusted to guarantee
the maximal signal from the sample. With surface muon beamlines (e.g. GPS at
PSI), the energy of the muons is usually 5 - 40 MeV/c. This is enough to
penetrate a thin layer of mylar foil or capton tape which is used to mount the
sample to the sample holder. However, this energy is not enough to penetrate
the thick layers of metal that make up a pressure cell. For this, decay muon
beamlines are needed, which can provide muons in the energy range of 40 - 125
MeV/c. For the surface muon case, alignment is simple as there is essentially
only the sample signal to detect. When the pressure cell is introduced, it
gets slightly more difficult as there is now another material in very close
proximity to the sample that will give a significant signal. This is not an
issue when the sample gives a very strong response, as the sample signal will
be much greater than the pressure cell signal. Consequently, determining the
optimal muon momentum for alignment is not very difficult. However, if the
sample response is weak, then it may be difficult to distinguish the sample
contributions from the pressure cell (as mentioned in Ref [8]). Additionally,
there are different types of pressure cells (MP35N, CuBe and a combination of
the two) with different signals. Further, each individual pressure cell made
from same material are in fact also slightly different. As a result, there are
several cases where it is incredibly difficult to extract the sample response
from the total signal.

One way around this is to first employ an indirect alignment using another
sample with a strong response (but same pressure cell), and then using this
alignment result also for the sample of interest (with weak signal). However,
this can be a problem as slight differences in the sample densities can
drastically change the stopping fraction of muons in the sample (a higher
density means less muons pass through the sample). This then translates into
issues in the fitting process as you need to have a correct estimate for the
signal fraction from the sample. Consequently, if your assumed stopping
fraction is incorrect, the obtain fitting results do not completely reflect
the intrinsic physical properties and behaviour of the sample.

In this paper, we present an efficient and user friendly simulation method
based using TRIM [9] combined with the pySRIM python module [1]. In such
approach, it is possible to simulate, with higher level of accuracy, the
number of muons that will be stopping in the pressure cell and the sample,
respectively, for any given muon momentum input. This is useful as a tool for
user to employ both before and after experiments. This is because it will both
give the user an idea on feasibility of their proposed experiment (in the
respect of how easy it will be to see their sample signal) and also for
supporting the data analysis when fitting the percentages of the contributions
from each fit function (pressure cell and sample).

##  2 Basis of Simulation

###  2.1 Beam setup

The GUI software is a python package that determines the stopping fractions in
a beamline setup. The software runs TRIM calculations in the background and
presents the results as early interpretable figures and stopping percentages.

There are three key aspects that make this software more accurate for
modelling stopping fractions compared to running the standard TRIM
simulations:

  1. 1.

The software utilises a Gaussian distribution of momentums/energies.

  2. 2.

The input beam can be collimated to any area.

  3. 3.

There is a random small angular divergence on each muon as each muon will not
be completely parallel to the beam.

Setting up the energy distribution in the beam first requires an input for the
percentage momentum acceptance. This represents the full width half maximum of
the momentum distribution. At the $\mu$E1 beam line, the momentum acceptance
is 3$\%$ [8], meaning a mean momentum of 95 MeV/c will give a full width half
maximum of 2.85 MeV/c. Thereafter a momentum range of 3$\sigma$ is selected on
both sides of the mean momentum. This is the initial input to the software
that then samples 20 discrete momenta within the $\pm$3$\sigma$ range. The
number of runs at each of these momenta is then weighted so that the sum of
these values gives a result as close to input number of runs, whilst also
following a Gaussian distribution, as seen in Fig. 1. The scaling value is
equal to the input number of runs divided by the sum of the selected
$y$-values on the Gaussian. This means that the total number of simulations
will either be less or more than the given total number of simulations.
However, it is never by such a large value that it becomes an issue. For 1000
runs (where one run is one muon simulated) at a momentum of 95 MeV and 3$\%$
momentum acceptance, 998 runs were completed, and for 10000 runs at the same
energy 1006 runs were completed. The momentum values then need to be converted
into their energy value, as TRIM only takes energy values as the input. We
then create a list of corresponding energy values using the following
equation:

$E=\sqrt{\left(pc\right)^{2}+\left(m_{0}c^{2}\right)^{2}}-m_{0}c^{2}$ (1)

where $p$ is the momentum, $c$ is the speed of light, and $m_{0}$ is the rest
mass of the muons $\left(m_{0}=105.658~{}MeV/c^{2}\right)$.

Figure 1: A graph to show how the number of runs at each energy is decided.
The height at each momentum is weighted such that it fits the Gaussian shape,
and when summed with the other heights it equals a value very close to the
number of simulations input into the software.

Each muons initial conditions (energy, position and angular divergence) is
written into the trim.dat file, which is provided when the SRIM software is
installed. The initial positions are decided according to the collimation at
GPD, 4x10mm2 [8]. Therefore the $y$ and $z$ positions are randomly sampled
according to a Gaussian distribution with full width half maximums of 2mm and
5mm respectively, to match beam size measurements taken at GPD. A small
angular divergence is added in also (-5∘ to 5∘ along all three axes) using
numpy.random.uniform in python. This in total created more accurate
representation of the beam setup at GPD.

###  2.2 Modelling the GPD setup

TRIM only allows to form layers of materials, and there is no possibility to
introduce either overlapping layers, or materials that are completely
surrounded by other materials. In other words, we had to model the GPD setup
as accurately as possible using only stacked individual layers of different
thickness and density. In the current setup for the simulations we primarily
target only the Variox-Heliox cryostat used at GPD (but different setups are
planned to be implemented later). The Variox is the first layer the muons see,
consisting off 0.1 mm Mylar, 0.2 mm Aluminium Nitrogen shield and 0.8 mm
Copper. Further, the Heliox consists of 2 mm thick Brass.

After penetrating the cryostat layers, the muon will reach the pressure cell
and finally the sample. At GPD, 3 types of pressure cells are used, \small{1}⃝
made entirely of MP35N alloy (Nickel - 35 , Cobalt - 35, Cromium - 20 ,
Molybdenum - 10, with a density of 8.43 g/cm3), \small{2}⃝ made entirely of
CuBe alloy (Berylium - 3, Copper - 97, with a density of 8.25 g/cm3) and
finally \small{3}⃝ a double-wall cell consisting of an inner chamber of CuBe
alloy and the outer part of the pressure cell made of MP35N alloy. All three
cells have an outer diameter 24 mm along with a sample space with diameter 6
mm and a height 12 mm. For the double-wall cell \small{3}⃝, the inner CuBe
chamber has an outer diameter 12 mm [10]. All cells are in the simulation
represented as slabs of: ’pressure cell - sample - pressure cell’.

###  2.3 Obtained Results

From the simulations, the $(x,y,z)$ position of every muon is collected. The
only exception to this is when the energy of the muon is so high that it
passes straight through the entire setup (all layers) and exits on the other
side. In this case, the muons stopping position is not saved (as it is out of
the range of calculation). By collecting the coordinates of each muon, the
stopping fraction for the number of muons inside the sample can be calculated.
The $(x,z)$ plane cuts across the circular area of the sample/pressure cell,
and the $y$ direction describes the height of the sample/pressure cell.
Choosing the origin to be the sample center, the absolute value of the $x$ and
$z$ coordinates of the moun is extracted, and evaluate if the square root of
the sum of the square of such coordinate ($\sqrt{x^{2}+z^{2}}$) is less than
the radius of the sample area (3 mm). If this condition is met and the $y$
coordinate is within -6 and 6 (for a 12 mm height with $z=0$ as the middle oft
the sample) then the muon is counted as being inside the sample. This number
divided by the total number of simulations will finally give the stopping
fraction.

There are a few assumptions made for this simulation. As mentioned earlier, we
can only simulate the GPD setup as ”slabs” of material. Therefore, the sample
height is not introduced in the simulations, and only comes in when working
out the sample stopping fraction. This means we have taken small cylindrical
region of the sample slab as the ”sample region”, and all around this region
as counting as the sample. This is a small issue due to most samples being put
into the beam having a lower density than the pressure cell, meaning that the
muons can move more freely in $zy$ plane in the sample slab than they would in
the real setup at GPD. However, it is estimated that the extra distance
travelled by the muons in this region will only have minuscule effect on the
accuracy of the simulations.

To ’field test’ the software We have also compared the simulations to the
alignment results on our previous GPD measurement of the CrCl3 material, which
used the same setup as the simulations. The stopping fractions were calculated
from the asymmetry of the sample in the fit, by taking the total possible
asymmetry and dividing this by the fit asymmetry. This percentage is taken to
be the stopping fraction of muons in the pressure cell. Here we see that the
simulations match the momentum the experiment was aligned on, at 100MeV/c. The
absolute values of the simulations do not match the experimental results
completely. This could be due to incorrect initial asymmetries in the sample
stopping fraction calculations, or an incorrect spread on the collimation.
However, as long as the shape (trend) of the results match the simulations it
is clearly reliable enough for alignment and background extraction.

Figure 2: The alignment measurements for CrCl3 plotted alongside the
simulation results. The sample was aligned using 100MeV/c, which the
simulation has shown to be the best momentum to use for the highest stopping
percentage.

##  3 Software GUI

The interface for running this software (see Fig. 3) has been set up to be
user friendly. The goal is to both aid in the alignment process as well as to
support data analysis via the calculations of the fractional contributions to
the signal. The steps for using the graphical user interface (GUI) are the
following (also noted in Fig. 3):

  1. 1.

Enter the desired number of energies you want to simulate in _Number of
Momenta_ , and click _Refresh Number of Momenta_. This should create a number
of text boxes under _Momenta (MeV/c)_ , where you fill in the momenta of
interest.

  2. 2.

Enter the number of elements in your sample in _Number of elements_ and click
_Refresh Number of Elements_. You should then fill out the element (Note that
this is the chemical symbol with the correct capitalised letters) and its
associated stoichiomtery.

  3. 3.

Fill out the sample density, the number of simulations, the full width half
maximum and also the directory that the SRIM application is on your machine.

  4. 4.

Finally, select which pressure cell you would like to use for the simulations
and then hit _Run_. Depending on the number of simulations and momenta you
have chosen, this could take different amount of time.

  5. 5.

When simulation is done, you can plot the data using _Plot Scatter_ or _Plot
Histogram_ , or you can print your results so that you can obtain the exact
positions and also the stopping fractions.

Figure 3: An example of how the Graphical User Interface (GUI) has been used
to calculate the stopping percentages of muons in a sample of TiSe2 at
momentums of 85 MeV/c, 95 MeV/c, 100 MeV/c, 105 MeV/c and 110 MeV/c, using the
completely MP35N pressure cell. The numbering relates to the list in Section
3.

As can be seen from the simulations (graphs in Fig. 3), we can obtain the
crucial information about all the muons in the simulations, both collectively
and as individual muons. We can also see that for higher momenta the spread of
muons parallel to the muon momentum direction is larger than the lower
momentum muons. This is because when a higher momentum is selected, these
simulations will contain a larger spread of momenta, hence, the stopping range
will be larger as well. Additionally, the spread perpendicular to the muon
momentum is also increased due to larger scattering effect. Such scattering
also affects the simulation time, i.e. as the input momentum is increased, the
time is takes to complete the given number of simulations increases also due
to the greater scattering. At the point when the energy is great enough that
some of the muons will pass through the entire pressure cell, the simulation
time reaches a maximum. As the momentum is increased beyond this point, the
whole process will take less and less time until all the muons pass through
the pressure cell (and then there is nothing more to simulate).

##  4 Conclusion and Outlook

We have developed a user friendly tool that allows for more accurate
simulations of the stopping fractions of muons for high-pressure experiments.
The options include all the different types of pressure cells provided by the
GPD beam line at the Paul Scherrer Institute (PSI) at any given muon momentum
using the Variox-Heliox cryostat. The Gaussian distribution of energies in the
simulated beam better represents the $\mu$E1 beam line when compared to the
single energies of the standard TRIM simulation. The GUI that has been made
make it easy for any user to find the stopping percentage of muons in their
sample prior to their experiment. This will aid in the alignment process as
well as allow the users to extract the background from the collected data. By
making this tool available to users it can both help writing their proposals
via a feasibility check of proposed experiment depending on how strong their
sample response will be. Further, it can also help to help the fitting process
via better estimate of signal fractions. We think that this will help solving
a key issue faced in high-pressure research in the $\mu^{+}$SR community, and
therefore push the field forward.

In the future, we would like to further improve this toolbox by extending
pySRIM beyond the ability to only make slabs of material. This will clearly
lead to an even more accurate representation of the beam setup at GPD. We
would also like to implement the possibility to select different sample
environments. This includes the different cryostats available at GPD, as well
as the new sample holders with different dimensions and thicknesses.
Additionally, extending this to surface muon beamlines would also be useful.
As stated previously, these beamlines have less of an issue with extracting
the sample response. However, this can still be a useful tool as you could
gain accurate sopping fractions for muons in layers of material, for example.
Furthermore, other simulation softwares, such as Geant4, would be useful to
implement alongside this in a similar fashion. This way we can get different
perspectives, and perhaps more accurate results.

##  5 Acknowledgments

This research is funded by the Swedish Foundation for Strategic Research (SSF)
within the Swedish national graduate school in neutron scattering (SwedNess),
as well as the Swedish Research Council VR (Dnr. 2021-06157), and the Carl
Tryggers Foundation for Scientific Research (CTS-18:272). Y.S is funded by the
Swedish Research Council (VR) through a Starting Grant (Dnr. 2017-05078). Y.S.
acknowledges funding from the Area of Advance-Material Sciences from Chalmers
University of Technology. G.S. acknowledges funding from the Chalmers X-Ray
and Neutron Initiatives (CHANS) grant and European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 884104 (PSI-FELLOW-III-3i). We thank Matthias Elender and Andreas
Suter for valuable discussions. We also thank Ola Kenji Forslun, Yuqing Ge,
Gaia De Berardino and Konstantinos Papadopoulos for the CrCl3 alignment data.

## References

  * [1] C. Ostrouchov, J.F. Ziegler, J.P. Biersack, and M.D. Ziegler.  pysrim-python tool for accessing SRIM package.
  * [2] O.K. Forslund, D. Andreica, Ya. Sassa, H. Nozaki, I. Umegaki, E. Nocerino, V. Jonsson, O. Tjernberg, Z. Guguchia, Z. Shermadini, R. Khasanov, M. Isobe, H. Takagi, Y. Ueda, J. Sugiyama, and M. Månsson.  Scientific Reports, 9(1), 2019.
  * [3] J.Sugiyama, D. Andreica, O.K. Forslund, E. Nocerino, N. Matsubara, Y. Sassa, Z. Guguchia, R. Khasanov, F.L. Pratt, H. Nakamura, and M. Månsson.  Physical Review B, 101(17), 2020.
  * [4] J. Sugiyama, W. Higemoto, D. Andreica, O.K. Forslund, E. Nocerino, M. Månsson, Y. Sassa, R Gupta, R. Khasanov, H. Ohta, and H. Nakamura.  Physical Review B, 103(10), 2021.
  * [5] V. Grinenko, D. Das, R. Gupta, B. Zinkl, N. Kikugawa, Y. Maeno, C.W. Hicks, H.H Klauss, M. Sigrist, and R. Khasanov.  Unsplit superconducting and time reversal symmetry breaking transitions in sr2ruo4 under hydrostatic pressure and disorder.  Nature Communications, 12(1), jun 2021.
  * [6] M. Majumder, R. Gupta, H. Luetkens, R. Khasanov, O. Stockert, P. Gegenwart, and V. Fritsch.  Spin-liquid signatures in the quantum critical regime of pressurized CePdAl.  Physical Review B, 105(18), may 2022.
  * [7] E. Kermarrec, J. Gaudet, K. Fritsch, R. Khasanov, Z. Guguchia, C. Ritter, K. A. Ross, H. A. Dabkowska, and B. D. Gaulin.  Ground state selection under pressure in the quantum pyrochlore magnet yb2ti2o7.  Nature Communications, 8(1), mar 2017.
  * [8] R. Khasanov, Z. Guguchia, A. Maisuradze, D. Andreica, M. Elender, A. Raselli, Z. Shermadini, T. Goko, F. Knecht, E. Morenzoni, and A. Amato.  High pressure research using muons at the paul scherrer institute.  High Pressure Research, 36(2):140–166, 2016.
  * [9] Ziegler JF, Biersack JP, and Ziegler MD.  Srim – the stopping and range of ions in matter.
  * [10] Z. Shermadini, R. Khasanov, M. Elender, G. Simutis, Z. Guguchia, K.V. Kamenev, and A. Amato.  A low-background piston–cylinder-type hybrid high pressure cell for muon-spin rotation/relaxation experiments.  High Pressure Research, 37(4):449–464, 2017.