**Principle.** The principle of muPS is similar to that of GPS. The muPS uses, instead of satellites, three or more detectors to form the reference coordinate (Fig. 1). The receiver detector defines the relative position within this coordinate by using the following relationship:

*L**i* 2 = (*x**i* *- x*p)2 + (*y**i* *- y*p)2 + (*z**i* *- z*p)2, (1)

where *L**i* is the geometrical distance between the *i*th reference detector located at (*x**i*, *y**i*, *z**i*) and the receiver detector located at (*x*p, *y*p, *z*p). Since the muon's traveling speed is approximately the speed of light in a vacuum (*c*0)19, the distance (*L**i*) between the reference detector and the receiver detector is determined by the muon's time of flight (TOF) between them. If the clocks attached to the reference detectors (reference clocks) and those attached to the receiver detectors (receiver clocks) are synchronized, this TOF is defined by the difference of the muon's arrival time between the reference and receiver detectors. Since *L**i*, *x**i*, *y**i*, and *z**i* are known, Eq. (1) has three variables; hence a ternary operation of Eq. (1) with three independent muon tracks will suffice. However, when conducting navigation with wireless muometric navigation system (MuWNS), the time synchronization between the reference clocks and receiver clocks is not automatically guaranteed. Actually, the stand-alone clock attached to the receiver detector is more likely to have a time offset error (Δ*t*) due to the frequency-drift of the crystal oscillators (OCXOs). Therefore, for MuWNS, Eq. (1) will have one more variable such that:

*L**i* 2 = (*x**i* *- x*p)2 + (*y**i* *- y*p)2 + (*z**i* *- z*p)2+*s*2, (2)

where *s* (= *c*0Δ*t*) is the pseudo-length that comes from the time offset at the receiver detector. Consequently, one variable was added to Eq. (1) and thus, at least four reference detectors will always be required for MuWNS. The positioning scheme shown in Eq. (2) holds only when we can assume that the time offset (Δ*t*) stays at the same constant rate or nearly constant rate within the time range required for positioning. The errors in determining *L**i* are defined as Δ*L**i* in this work. For simplicity, it was also assumed that the reference detectors were constantly connected to the the GPS antenna and getting an uninterrupted signal.

**Required accuracy for MuWNS.** In order to evaluate what the most practical positioning precision will be for each given situation, referencing the accuracy of commonly used portable GPS receivers is useful since the studies reporting the accuracy of these kinds of GPS receivers are extensive. For example, von Watzdorf and Michahelles34 reported that the average accuracy of location information ranged between 108 and 655 m with GPS enabled Apple iPhones, iPods, and iPads. Zandbergen35 reported that the average horizontal position error of an iPhone was around 10 m, and Garnett and Stewart36 found better precision for positioning with an iPhone 4S which was around 6.5 m. However, these varieties of the positioning accuracy depend on the environment (such as how dense the buildings and trees, or topography are) since large structures surrounding the GPS receiver can cause multipath errors. In this work, we define the following three positioning precision categories for practical muometric navigation: Very Precise Level (VPL), a precision better than 1 m (Δ*L**i* ≤ 1 m); Precise Level (PL), a precision better than 3 m (Δ*L**i* ≤ 3 m); and Usable Level (UL) a precision better than 10 m (Δ*L**i* ≤ 10 m). Positioning accuracy worse than UL will be defined as an unpractical level and will not be discussed in this work.

**Design of MuWNS.** Figure 2 shows the block diagram of the current MuWNS. The time of the MuWNS setup is defined by the GMC, and every MuWNS detector records absolute time information. The GMC is connected to the GPS antenna when GPS signals are available. When GPS signals are not available (for some reason such as multipath reception, GPS jamming and spoofing, system failures, diving underwater, moving underground, etc.), the GMC will automatically switch from the GPS mode to the holdover mode. The muPS signals (e.g. photomultiplier tube (PMT) outputs) are discriminated with the constant fraction discriminator (CFD) and the time difference from the GMC signal is measured with the TDC (Fig. 2A). When GPS signals are not available, the event-counter-coincidence (ECC) unit starts to count the GMC signals at a rate of 10 MHz. When the GMC and muPS coincidence signals are generated, the number of GMC counts (*n*) and the timing information (*t*) are transferred to a microcomputer (e.g., Raspberry Pi) respectively from ECC and TDC. The transferred *n* and *t* values are shared between MuWNS detectors via Wi-Fi or an acoustic modem to calculate the TOF. As a consequence, each MuWNS detector records the muon arrival time in a sharable form: *n*T + *t*, where T is the GMC/holdover signal period, and the "time zero" is defined as the moment when the GPS signal is lost (Fig. 2B); hence the TOF between the reference detector and the receiver detector can be derived by comparing reference and receiver clocks: |[*n*T + *t*]REFERENCE-[*n*T + *t*]RECEIVER|. The TDC is designed to reset start signals if it doesn't receive a stop signal before the next start signal, and the TDC time range must be sufficiently longer than T (in this case 100 ns). According to Tanaka (2020)19, since the jitter level for *t* is less than 200 ps, fluctuations in T will be a major factor to degrade the current TOF measurement quality. The errors in T will be described in the following subsections.

**Grandmaster clock stability.** The Grandmaster clock (GMC) stability is the key component to define the accuracy of MuWNS. Experimental results of the current GMC/holdover are described here. The timing jitter measured with the current GMC is shown in Fig. 3. These time-dependent variation plots were drawn in the following way. The OCXO used in the current work (Trimble OCXO) was equipped inside the GPS grandmaster clock (Trimble Thunderbolt PTP GM200). In this study, the Trimble OCXO was initially synchronized with the GPS signal by connecting a GPS antenna (See Method section). The OCXO was then subsequently disconnected to simulate indoor/underground/underwater navigation. This process was repeated 12 times to investigate the drift pattern. The results showed that the drift speed varied in each run, and both positive and negative drifts were measured. If the OCXO operation time lasted less than 1 hour, the drift was roughly linear as a function of time, but during longer operations, the behavior was unpredictable. Figure 3 shows the averaged time offset (the deviation between the GPS and holdover mode averaged over the entire run: <Δ*t* (t)>) and the standard deviation (SD) of Δ*t* (σ(*t*)) as a function of time. Since both the drift speed and the drift direction were unpredictable, σ(*t*) will remain as a non-zero factor. In MuWNS operations, Δ*t* can be derived by solving Eq. (2), but non-zero σ will lead a major positioning error.

**Time synchronization between reference clocks and underground/underwater clocks.** In order to evaluate the time synchronization capability in Eq. (2) between the reference and receiver detectors with the currently designed MuWNS, the relationship between the GMC stability and the muon tracking frequency were studied. If the muon tracking frequency was too low, the uncertainties in the drift level and direction would increase, and σ(*t*) would not be negligible, and as a consequence, quartic equations in Eq. (2) cannot be approximated as simultaneous equations anymore.

The muon's tracking frequency (*N*− 1) can be derived by integrating the open-sky zenith angular muon's energy spectrum *I* (*E*, *θ*)37,38,39 over the energy range between *E*c and infinity:

*N*(*θ*)=\({\int }_{{E}_{c}}^{\infty }I(E,\theta )dE\), (3)

where *E*, *θ*, *φ* are respectively the muon energy, elevation angle, and azimuth angle. *E*c is the muon cutoff energy that depends on azimuthal variations of matter thickness located above the detector. For a derivation of *E*c, the muon's constant slowing down approximation (CSDA)40 can be used. In order to keep the relativistic energy of muons in a state sufficiently higher than their invariant mass, *E*c in Eq. (3) was set to be 1 GeV higher than the cutoff energy directly derived from the CSDA muon range. *N* in Eq. (3) is also a function of the solid angle (*Ω*) formed by the reference and receiver detectors. If the size of the detectors (*S*) is given, *Ω *depends on the distance between the reference and receiver detectors.

**Required spatial density of the reference detector network and resultant navigation performance.** In order to evaluate spatial density of reference detectors required for achieving a given level of navigation performance of the currently designed MuWNS, the positioning accuracy and the time required for attaining this accuracy were studied. Generally, in the scheme of muPS, recording more muon tracks will improve the positioning accuracy more since the muon's arrival positions to the detector would be better averaged over the detection area, thus more muon events provide better positioning precision19. In the MuWNS scheme, the positioning errors coming from SD of Δ*t* have to be considered in addition to the positioning errors coming from this detector size effect. Since σ increases as time goes by, the muon's tracking frequency and the positioning accuracy is tradeoff. In the following discussion, we assume the detector size is negligible in comparison to Δ*L**i*; hence Δ*L**i* = *c*0σ(*t*).

Figure 4 shows *c*0σ(*t*) as a function of the maximum distance between the reference and receiver detectors, where the maximum distance was defined as the distance (i.e., *SΩ*) that enables each pair of the reference and receiver detectors to record at least one muon track within the given time so that *c*0σ(*t*) is less than the given accuracy level. In other words, if the distance between the reference and the receiver detectors exceeds this maximum distance, the given positioning accuracy will not be attainable due to the limited flux of cosmic-ray muons. Therefore, the maximum distance defines the required density of the receiver detectors for the required positioning precision. In this calculation, the effective active area of the reference and receiver detectors were both assumed to be 1 m2, and a typical basement floor without surface structures above it was assumed; namely, the assumed characteristics were: 5-m thick soil (SiO2) with an average density of 2.0 gcm− 3 existing between the reference and the receiver detectors. As the muon's arrival angles approach the horizontal angle, the flux is reduced. As a consequence, the distance between the reference and receiver detectors has to be shortened to attain larger *Ω*. In this case, it is necessary to place reference detectors more densely.

Although the OCXO's frequency drift depends on a complicated combination of the environmental factors, it was expected that the resultant σ(*t*) would be suppressed if we utilize multiple independent OCXOs and take an average of σ(*t*) between them. Figure 5 compares σ(*t*) averaged over single, double and triple OCXO outputs. For example, if we use three OCXOs, the time we are allowed to keep the low drift level to attain VPL will be extended from less than 9 seconds to less than 28 seconds. As a result, a less dense reference detector array will be required for a multiple OCXO system.