Optimization of muonium yield in perforated silica aerogel

A muonium consists of a positive muon associated with an orbital electron, and the spontaneous conversion to antimuonium serves as a clear indication of new physics beyond the Standard Model in particle physics.One of the most important aspects in muonium-to-antimuonium conversion experiment (MACE) is to increase the muonium yield in vacuum to challenge the latest limit obtained in 1999. This study focuses on a simulation of the muonium formation and diffusion in the perforated silica aerogel. The independent simulation results can be well validated by experimental data. By optimizing the target geometry, we find a maximum muonium emission efficiency of $7.92(2)\%$ and a maximum vacuum yield of $1.134(2)\%$ with a typical surface muon beam, indicating a 2.6 times and a 2.1 times enhancement, respectively. Our results will pave the way for muonium experiments.


I. INTRODUCTION
Muonium, as a pure leptonic bound state (M = µ + e − ), is a ideal subject in precision measurement of the Standard Model and a sensitive probe to new physics beyond SM.
Improving the yield of vacuum muonium will make significant contributions in muonium spectroscopy [1][2][3][4][5], precise measurement of muon g−2 and electric dipole moment (EDM) [6][7][8], muonium antimatter gravity experiment (MAGE) [9][10][11], and the search for muoniumto-antimuonium conversion [12][13][14][15].The conversion process has not been studied by any experiment since 1999 and the latest limit on the conversion probability is P MM < 8.3 × 10 −11 (90% C.L.) [16].Recently, we proposed a new muonium-to-antimuonium conversion experiment, MACE, aiming at improving the sensitivity to the spontaneously charged lepton flavor violating conversion process (M → M, or µ + e − → µ − e + ) by more than two orders of magnitude [17,18].MACE is designed to search for the rare process by distinguishing charge and kinematic properties of the decay products between muonium and antimuonium, where only vacuum muonium conversion events are of great interest.These muonium-related experiments are benefited strongly from the increase of the muonium yield in vacuum.
How can we produce vacuum muonium?A muon beam is injected into a specific material and a muon may spontaneously capture an electron to compose a muonium atom.After that, some of produced muonium atoms would eventually escape from the material and emit into vacuum.The target material is chosen to be porous and inactive such as silica powder or silica aerogel.As a matter of fact, the silica powder was used in the latest experiment searching for muonium-antimuonium conversion [16].However, the silica powder is fragile and has relatively low vacuum muonium yield.The silica aerogel shares the similar composition but takes a rich porous meso-structure, suitable for the muonium target.Early in 1992, the muonium emission efficiency in a silica aerogel was measured for the first time and the result was comparable to that of the silica powder, even though a muon beam with higher momentum was used (which is likely to be adverse to muonium emissions into vacuum as discussed in section IV) [19].After 30 years, a new measurement was carried out and the result demonstrated statistically that up to 3 muonium atoms could emit into vacuum out of 1000 stopped muons [20].Keep in mind that most of the muonium atoms formed in the material are still trapped [20][21][22][23].There remains a plenty of room for improvement in the muonium emission into vacuum and a new approach is expected.G. A. Beer et al. proposed a laser-ablation approach aimed at increasing the muonium emission efficiency.In this approach, an aerogel target is ablated with a femtosecond pulsed laser so that the target surface at the beam downstream side is perforated (Fig. 1) [24].Laser ablates nearly cylindrical holes with diameters ranging from tens to hundreds of micrometers and depths of a few millimeters.The holes are arranged in an equilateral triangular lattice, perpendicular to the surface of the target with a spacing of tens to hundreds of micrometers.
The perforated structure significantly enhances the diffusion of muonium and dramatically increases its vacuum yield.Multiple ablation structure parameters were tested in TRIUMF and an emission efficiency of up to 2% was achieved, indicating an order of magnitude improvement [24].The authors also attempted to reveal the emission mechanism of muonium by proposing and analyzing the opening fraction, a characteristic parameter, to figure out the dependency of muonium emission efficiency.However, the latest results roughly showed a positive correlation between the opening fraction and yield while more complex dependencies remained blank.The insufficient number of samples also led to an incomplete understanding of the underlying mechanism.In addition to the experiments conducted at TRIUMF, measurements on muonium emission efficiency in perforated aerogel were also conducted in PSI [25].In the PSI experiment, a slow muon beam (p µ + = 11-13 MeV/c) entered the front surface of the aerogel sample, and muonium atoms were emitted into vacuum upstream.However, the emission efficiency depends on the emission model parameters and further discussion on the enhancement of muonium emission by perforated structures is missing.We also note that an independent simulation study was carried out to model the muonium diffusion in silica aerogel and was applied to a novel multi-layer target design [26].Their study modelled the transport of muonium as a diffusion process with a thermal momentum distribution.However, they incorporated the effect from enlarged diffusion constants while detailed considerations of perforated structures remained blank.In this study, we will model muonium diffusion as a three-dimensional off-lattice thermal random walk and propose a new muonium tracking algorithm that comprehensively considers the perforated structure.We will simulate the muonium diffusion inside perforated silica aerogel with the measured mean free path rather than fitting the diffusion constant so that we will carry out the first explicit dependency of the muonium yield on the perforation parameters.Furthermore, the simulation method will be validated and applied to a muonium yield optimization.We will attempt to analyze and understand the dependencies and the underlying mechanisms involved in the muonium emission.This study will further guide the development and optimization of the muonium target, and pave the way for the design and optimization of MACE.Our method will not only benefit the precision measurement and the search for new physics in the muon sector, and the development on cutting-edge muon spin relaxation/rotation/resonance (µSR) techniques [27][28][29][30][31], but also be meaningful to the positronium physics and technologies, including the positronium production and emission [32][33][34][35], positronium spectroscopy and annihilation [36][37][38][39], antihydrogen gravity experiment (AEḡIS) [33,40,41], and positron emission tomography (PET) [42].
This article is organized as follows: In section II, the theoretical formalism and the simulation method are described.In section III, we present the validation of the simulation.
We then apply the new algorithm to optimize the vacuum muonium yield in section IV.The findings of this work are finally summarized in section V.

II. MUONIUM DIFFUSION AND ITS SIMULATION
Silica aerogel is an inactive porous material with a rich meso-structure.It is a matter similar to common gels, but liquid component are replaced by gas, which makes its porosity extremely high (80%-99.8%)but density extremely low (O(mg/cm 3 )) [43][44][45].Microscope photos demonstrate that a silica aerogel has a complex network structure, with a characteristic scale of tens to hundreds of nanometers.The silica network structure allows molecules or atoms diffuse inside the inter-granular void [46].
If the aerogel is placed in vacuum for a sufficient period, air molecules inside would diffuse and come out of the material.The aerogel would become a "vacuumgel" and muonium can move freely inside the inter-granular void until colliding with silica.After a collision, muonium will mostly be inelastically scattered instead of being captured.This proposition is supported by density functional theory (DFT) calculations of hydrogen atom dynamics inside silica.Based on DFT results, we can infer that a majority of muonium atoms can disassociate from silica and remain its state neutral [47].Independent experimental results on absorption of hydrogen by silica aerogel also demonstrate that significant absorption of hydrogen only takes place at a very low temperature, supporting this proposition to a certain extent [48].Therefore, muonium atoms experience a cycle of free motion -scattering -free motion process in a silica aerogel before decays and it can be modelled as a three-dimensional  Following this principle, a recursive equation of motion describing a muonium atom inside a silica aerogel can be written as where (t i , x i ), r i , v i and v i = ∥v i ∥ are the space-time coordinates, the free path, the velocity and the velocity magnitude at the i-th step, respectively.Both r i and v i are random variables sampled from the corresponding distribution at each step.Considering the randomly connected network structure in a silica aerogel, it is convinced that the free path follows an exponential distribution with a parameter λ where λ = ⟨r⟩ is the mean free path depending on the structure of a specific aerogel.The magnitude of λ is positively correlated with the size of meso-cavities inside the aerogel.
It comes to the velocity distribution now.Inside an aerogel, muonium atoms experience a series of inelastic scattering on silica structures, and the emission velocity distribution after scattering is related to the interaction between the muonium atoms and silica.[47,49,50].
Therefore, in principle, the exact distribution should be determined by experimental data or ab inito calculations.As an approximation, it is believed that the ensemble of muonium atoms and target molecules quickly reaches the thermal equilibrium [49].Therefore, the kinetic energy of muonium atoms should follow the form of the Boltzmann distribution, and the velocity distribution takes the form where m M is the mass of muonium, k B is the Boltzmann constant and T eff is the effective temperature.It is known that the room temperature will overestimate the muonium decay time.One of the interpretation is the underestimation of kinetic energy.We can always correct the effect by elevating the temperature [20,26].The experimental work of Antognini et al. found the temperature to be 400 K, by considering a simplified emission model with a given diffusion time ∆t diff = 200 ns [25].However, their emission model has correlations between ∆t diff and temperature.Therefore, the temperature depends on the choice of ∆t diff .
Another simulation study presented by Ce Zhang et al. estimated the temperature to be around 320 K (320 ± 2 K for flat target and 322 ± 2 K for perforated target) with datasets measured at TRIUMF.In the model, the diffusion time is naturally included and the temperature is determined by data fitting.In out study, we choose T eff = 322 K to perform the following simulations.Equation ( 1)-( 3), as an iterative scheme of the special off-lattice random walk, simulates the motion of muonium until it reaches the material surface or decays.
To perform accurate simulations in a perforated structure, it is also important to track muonium outside the material reliably.The commonly used method for navigating and tracking particles in a complex geometry is to calculate the geometric relations between particles and the voxelized geometry (e.g.Geant4 [51][52][53]).However, this method can hardly handle a perforated target with tens of thousands of micro holes, since neither the full geometry can be easily constructed nor an efficient navigation or tracking can be performed.Therefore, an alternative method is needed.Our solution is to use a combinatorial tracking method.The voxelized-geometry-based tracking method controls muonium motion outside the target bounding box while a Boolean-expression-based tracking method controls muonium motion inside.This combination can resolve the realizability and performance issues raised in the voxelized-geometry-based tracking method.
To implement the method, the first question is how to describe the target full geometry in the combinatorial tracking system.The solution involves a Boolean expression to describe the perforated structure and a bounding box This combinatorial tracking method combines the power of Geant4 and the Booleanexpression-based algorithm in the same procedure, allowing us to successively simulate a muon track and a muonium track.We are able to simulate the dependency of the muonium yield with regards to different muon beam conditions conveniently and precisely, so that the approach can be validated with experimental data.The large-scale simulation can be performed in a local supercomputer following the software deployment strategy in Ref. [54].

III. VALIDATION OF SIMULATION SCHEME
After a discussion of how to simulate the muonium diffusion in a perforated target, this section demonstrates a validation of the simulation method.The previous experimental studies, utilizing aerogel targets with many tiny blind holes arranged in an equilateral triangular lattice at the downstream surface, have accumulated vacuum muonium production data in different experimental setups.Two datasets can be used for the validation: one is the vacuum muonium yield data [23], and the other is the spacetime distribution of muonium decays [20,24].Simulations will be performed and corresponding datasets will be derived  by analyzing the simulation data.
The model described in section II contains two free parameters: a mean free path λ and a temperature T .While the temperature is chosen to be 322 K according to the previous fit [26], the mean free path is estimated by the empirical scaling formula λ = λ 0 (ρ 0 /ρ) 1.5 [20].The reference mean free path λ 0 and the reference target density ρ 0 are 226 nm and 29 mg/cm 3 , respectively.Since the target densities involved in this validation are close to the reference values, the scaling formula can be applied.The perforated structure is simplified as cylindrical blind holes arranged in an equilateral triangular pattern, with its axis perpendicular to the target surface, as shown in Fig. 1.This structure is represented completely by a Boolean expression and a target bounding box (Geant4 G4Box).The spacing, diameter, and depth of the holes are geometric parameters kept the same as the experiment.The remaining simulation settings are also kept the same as the corresponding benchmark experiment.
We first validate the spacetime distribution of muonium decay by comparing it with measurements in TRIUMF [20,24].In the simulation, a σ xy = 5 mm Gaussian profiled collimated surface muon beam carrying 22.9 MeV/c momentum with 1.5% spread (rms) is incident on the target.The target thickness is 6.9 mm and perforated by 270 µm diameter and 30 µm spacing holes, from which some of the muonium atoms emit into the surrounding vacuum and decay.These decays are divided into three different space regions ranging from z = 10 mm to 40 mm, as shown in the simulation results (Fig. 3).The simulated distribution and an exponential muon decay background are fitted to the experimental histogram.A room temperature (293 K) and an elevated temperature (322 K, according to Zhang's fit [26]) are used in the simulation.The result shows agreement between simulation and experiment with χ 2 (322 K) < χ 2 (293 K) for all cases, as we mentioned above.We believe the goodness of fit can be further improved by tuning temperature for each case exclusively.In this validation case, we deal with the decay spacetime distribution of muonium emission.Our primary goal is to simulate the vacuum yield of the perforated target.One might ask whether the small difference in the decay spacetime distribution will significantly affect the yield result.We will answer this question by a second round of validation on the vacuum muonium emission efficiency.
In the simulation, a collimated surface muon beam with 23 MeV/c momentum and 2% hardly covered in the simulation.We tend to consider the deviations as the systematic errors, which will not affect the reliability of the simulation.In conclusion, the straightforward simulations now can describe the transport of muonium well in complex aerogel geometry and provide reasonable yield estimations.Precision measurement experiments might need calibrated simulations in details.

IV. SIMULATION-GUIDED YIELD OPTIMIZATION
With the simulation method validated, we can utilize it to investigate the performance of perforated targets of different ablation configurations to guide the optimization of the mainstream scheme in the related muonium experiment.In this section, we will study the muonium yield under various target geometries to find the dependencies of muonium yield on key parameters.The perforation structure is the same as Fig. 1, with a target size of  and 10% are considered.In a real experiment, the muon beam usually passes materials such as a beam degrader or a beam monitor [55].Therefore, we place an aluminum degrader at 5 mm in front of the target where the muon beam goes through and loses energy.A portion of muons will stop in the target and form muonium atoms.The degrader with a Ref. [25] (Aerogel-1) 12.5 3. We use five quantities, f stop µ + , f M , Y tot , R vac , and Y vac , to characterize the formation and emission of muonium.They represent the muonium formation fraction, the muon stopping fraction, the muonium total yield, the muonium emission efficiency, and the vacuum Ref. [25] (Aerogel-1) a 12.5 3. b Converted from R vac (Table I) by muonium yield, respectively.They are defined as Among them, N stop µ + is the number of muons stopped in target, N in µ + is the total number of incident muons, N tot M is the total number of produced muonium, and N vac M is the number of vacuum muonium.There are following notable relationships between them While f M is mostly determined by material properties and f stop µ + is related to both material properties and the beam condition, other three quantities, Y tot , R vac , and Y vac , will be affected by aerogel target geometry.On one hand, increasing the spacing between holes will increase the total amount of material, thereby increasing the rate of muons stopping in the material and producing muonium, but these muonium atoms are less likely to emit into vacuum due to the material's obstruction.On the other hand, increasing the diameter of the hole can improve the emission efficiency, but this will reduce the overall muonium yield.Therefore, it can be expected that there is a geometric parameter combination that leads to the highest vacuum yield when the material and the beam condition are fixed.The simulation results support this prospect.As shown in Fig. 6, there exists a maximum value for both vacuum yield and emission efficiency.As the spacing and diameter deviate from the optimal values, the vacuum yield and emission efficiency gradually drop and are significantly suppressed when the spacing or diameter is too tiny.
Based on the notable relation Y vac = R vac f M f stop µ + , the enhancement in muonium vacuum yield Y vac can be achieved by three different mechanisms.
• Enhancement by sufficient open surface area that enhances muonium diffusion towards vacuum, which contributes positively to R vac .This enhancement is suppressed when the diameter is too small or the spacing is too large, as shown in Fig. 6.In this case, although muonium atoms are formed near the target surface, there is little opening surface for them to easily diffuse into the vacuum.Meanwhile, they have a long distance to travel before leaving the material, during which most of them decay and cannot reach the target surface.
• Enhancement by favorable forming position.In this case, muoniums originates near the target surface and is favorable for their emission into vacuum, which contributes positively to R vac .This enhancement is suppressed when the spacing is tiny, as shown in Fig. 6, or the beam momentum spread is large, as shown in Table I and Table II.For tiny spacings, muonium tends to form near the bottom of the hole rather than near the target surface, since there is not enough material to support sufficient muonium formation near the surface.For large beam momentum spread, the distribution of muonium forming position tends to be uniform, with a large number of muoniums produced deep inside the target.Both of these two situations require a long diffusion path for muonium atoms to travel before they emit into vacuum.
• Enhancement by appropriate material removal, which contributes positively to Y tot .This mechanism takes effect when materials are not excessively removed (i.e. with large spacing or small diameter, as shown in Fig. 6).At this point, there are more muons stopped inside the target so the original total yield of muonium is higher, leading to a higher vacuum yield.
In conclusion, we prefer to have a high-quality muon beam with minimize momentum spread and an appropriate opening surface area, which can maximize the number of muonium atoms forming near the target surface and maximize the muonium emission into vacuum.This guides us to optimize the muonium target in the experiment MACE.

V. SUMMARY AND OUTLOOK
We have developed a technique to simulate the muonium production and tracking in perforated silica aerogel, and implemented it in Geant4 as an independent physical process.This simulation technique has been validated with the spacetime distribution in muonium decays and the muonium emission efficiency measured in TRIUMF.The validation has shown a good consistency between simulation and experiment.In addition, optimization of geometric parameters in the perforated aerogel has been carried out so that we could make good use of the simulation technique to facilitate muonium-related experiments such as MACE.We have considered a simplified simulation setup based on real experimental circumstances and incident beams under different conditions to investigate the performance of different target geometries.We have scanned the depth, spacing, and diameter of the holes within a certain range and, for the first time, have determined the explicit dependency of the muonium total yield, the muonium emission efficiency, and the vacuum muonium yield on the geometric parameters.We have identified the mechanisms for the enhancement of vacuum muonium yield in the experimental design and possible optimization directions are indicated.Nevertheless, there is still plenty of room for further improvement in the simulation techniques and muonium targets.On one hand, in the simulation technique, the effect of material removal caused by laser ablation on the muon transportation has not been considered in detail due to the independence of the current process implementation, and the meticulous combination of this muonium tracking process with the muonium conversion process or other muonium-related processes inside the material is still a technical topic needs to be studied.On the other hand, it is convinced that the combination of a low-momentumspread muon beam and an optimized target will lead to higher vacuum muonium yields, which deserves further investigation.We believe that this new simulation method can inspire novel target designs and help achieve better results.It could be expected that our simulation technique will make an essential contribution to the optimization and guidance of muonium target design and will play an important role in the related science and technology frontiers, which share similar hydrogen-like atom formation and diffusion physics.

FIG. 1 :
FIG. 1: The schematic diagram of muonium emission in the silica aerogel target.

FIG. 2 :
FIG.2:The muonium formation and tracking algorithm in the simulation.
(a Geant4 volume) to describe the target region where the Boolean expression should be enabled.The Boolean expression divides the target volume into two different sub-regions: a true region representing the material where muonium atoms follow the random walk, and a false region representing the vacuum where muonium atoms move freely.During each step of the random walk, the Boolean expression is calculated to determine the inclusion of the muonium in the material and the state of motion.If a muonium atom appears inside the material, a free path and a velocity are sampled and the spacetime coordinates are updated.Otherwise, only a free path is sampled and the state of motion is updated with the current velocity.The process will loop until any of the termination conditions, i.e., the muonium decays or escapes from the target bounding box, is triggered.Then the Boolean-expression-based tracking comes to an end and the tracking procedure is returned to the voxelized-geometry-based method (Geant4).The muonium formation and tracking algorithm in the simulation is shown in Fig.2.

(f) Perforated target, region 3 FIG. 3 :
FIG. 3: The spacetime distribution of muonium decay.The simulation histogram and an exponential muon decay background are fit to the experimental histogram.

FIG. 4 :
FIG.4:A comparison of muonium emission efficiency in simulation and experiment.

FIG. 6 :
FIG.6: Projection of muonium total yield, muonium emission efficiency, and vacuum muonium yield under different beam conditions and target geometries.

4 4 . 5 ±a
0.5 45 ± 5 105 ± 5 6.72 ± 0.05 +1.06 −0.76 c Only statistical errors are shown.b Includes muonium decays within 10 mm < z < 40 mm only.c Model-dependent assumptions of the temperature at 400 K and the diffusion time at 200 ns.differentthickness will result in various distributions for muon stopping positions.Thereby, the degrader can affect the emission of muonium into vacuum and eventually change the vacuum yield.It means that the degrader thickness should be optimized to achieve the maximum yield, as shown in simulation results Fig.5.From the results, we infer the optimal degrader thicknesses for beam momentum spread of 2.5%, 5%, and 10% are 430 µm, 410 µm, and 370 µm, respectively.

√ 3 / 2 (
s+ϕ) 2 by considering muons are half-stopped in flat target and taking into account the material volume fraction near the surface after ablation, and f = 0.655 [23].c Includes muonium decays within 10 mm < z < 40 mm only.d Model-dependent assumptions of the temperature at 400 K and the diffusion time at 200 ns.

TABLE I :
Simulation of maximum muonium emission efficiency and corresponding optimal spacing and diameter with different beam condition.

TABLE II :
Simulation of maximum vacuum muonium yield and corresponding optimal spacing and diameter with different beam condition.