# No Large Dependence of Planet Frequency on Galactocentric Distance

Naoki Koshimoto  Laboratory for Exoplanets and Stellar Astrophysics,
NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA  Department of
Astronomy, University of Maryland, College Park, MD 20742, USA  Department of
Earth and Space Science, Graduate School of Science, Osaka University,
Toyonaka, Osaka 560-0043, Japan  David P. Bennett  Laboratory for Exoplanets
and Stellar Astrophysics, NASA/Goddard Space Flight Center, Greenbelt, MD
20771, USA  Department of Astronomy, University of Maryland, College Park, MD
20742, USA  Daisuke Suzuki  Department of Earth and Space Science, Graduate
School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan  Ian A.
Bond  Institute of Information and Mathematical Sciences, Massey University,
Private Bag 102-904, North Shore Mail Centre, Auckland, New Zealand

###### Abstract

Gravitational microlensing is currently the only technique that helps study
the Galactic distribution of planets as a function of distance from the
Galactic center. The Galactic location of a lens system can be uniquely
determined only when at least two of the three quantities that determine the
mass–distance relations are measured. However, even if only one mass–distance
relation can be obtained, a large sample of microlensing events can be used to
statistically discuss the Galactic distribution of the lenses. In this study,
we extract the Galactic distribution of planetary systems from the
distribution of the lens-source proper motion, $\mu_{\rm rel}$, for a given
Einstein radius crossing time, $t_{\rm E}$, measured for the 28 planetary
events in the statistical sample by Suzuki et al. (2016). Because microlensing
is randomly caused by stars in our Galaxy, the observational distribution can
be predicted using a Galactic model. We incorporate the planet-hosting
probability, $P_{\rm host}\propto M_{\rm L}^{m}R_{\rm L}^{r}$, into a Galactic
model for random-selected stars, where $M_{\rm L}$ is the lens mass ($\sim$
host mass), and $R_{\rm L}$ is the Galactocentric distance. By comparing the
observed distribution with the model-predicted $\mu_{\rm rel}$ distribution
for a given $t_{\rm E}$ at various combinations of $(m,r)$, we obtain an
estimate $r=0.2\pm 0.4$ under a plausible uniform prior for $m$ of $0<m<2$.
This indicates that the dependence of the planet frequency on the
Galactocentric distance is not large, and suggests that the Galactic bulge
does have planets.

Exoplanets (498), Milky Way Galaxy (1054), Galactic bulge (2041),
Gravitational microlensing (672)

##  1 Introduction

Although more than 3000 planetary systems have been discovered to date, most
reside at distances $<1~{}{\rm kpc}$ from the
Sun111https://exoplanetarchive.ipac.caltech.edu/ (catalog NASA Exoplanet
Archive) (Akeson et al., 2013). In this regard, gravitational microlensing is
a unique technique because it is sensitive to planetary systems over a wide
range of distances in our Galaxy, from the Galactic disk to the Galactic bulge
(Gaudi, 2012). Microlensing is currently the only technique capable of
investigating the Galactic distribution of planets. However, no study
measuring the distribution has yet been reported, mainly owing to the
difficulty in distance measurement.

There are four physical quantities involved in each microlensing event: the
lens mass $M_{\rm L}$, distance to the lens $D_{\rm L}$, distance to the
source $D_{\rm S}$, and lens-source relative proper motion $\mu_{\rm
rel}$—–which is given by $\mu_{\rm rel}=|\mu_{\rm L}-\mu_{\rm S}|$, where
$\mu_{\rm L}$ is the lens proper motion, and $\mu_{\rm S}$ is the source
proper motion. For most planetary events, two parameters related to these
quantities can be commonly measured via light curve modeling, i.e., the
Einstein radius crossing time,

$\displaystyle t_{\rm E}=\frac{\theta_{\rm E}}{\mu_{\rm rel}},$ (1)

and lens-source relative proper motion, $\mu_{\rm rel}$. Here, $\theta_{\rm
E}$ is the angular Einstein radius given by

$\displaystyle\theta_{\rm E}=\sqrt{\kappa M_{\rm L}\pi_{\rm rel}}$ (2)

where $\kappa=8.144~{}{\rm mas}~{}M_{\odot}^{-1}$ and $\pi_{\rm rel}={\rm
1\,AU}(D_{\rm L}^{-1}-D_{\rm S}^{-1})$.

Measuring either the microlens parallax or lens brightness, besides $t_{\rm
E}$ and $\mu_{\rm rel}$, is required to determine $M_{\rm L}$ and $D_{\rm L}$,
even assuming that the source star is located in the bulge (i.e., $D_{\rm
S}\sim 8~{}$kpc). This requirement adds to the complexity of statistical
studies focusing on determining the Galactic distribution of planets based on
measured distances to the lens. The microlens parallax is measured only when
the lens is relatively close to the Sun (Bennett et al., 2010). Because of the
bias toward closer lenses and the vulnerability to systematic errors (Penny et
al., 2016; Koshimoto & Bennett, 2020), much effort is required to prepare a
clean statistical sample of a sufficient number of microlens parallax
measurements.

Lens brightness measurements seem to be more robust in terms of susceptibility
to systematic errors. However, to resolve the lens from the source, high-
angular-resolution follow-up observations with adaptive optics or observations
by the Hubble Space Telescope (HST) are required several years after the
event. Thus, a certain amount of time is necessary to obtain a statistical
sample of the lens brightness measurements. A systematic follow-up program for
measuring the lens brightness for past planetary events is ongoing using the
Keck telescope and HST (Bennett, 2018; Bhattacharya et al., 2018), and a
statistical sample of past planetary events with measurements of $t_{\rm E}$,
$\mu_{\rm rel}$, and the lens brightness will be available in the near future.

Meanwhile, this study aims to extract information about the Galactic
distribution of planets from a set of measurements of the Einstein radius
crossing time, $t_{\rm E}$, and the lens-source relative proper motion
$\mu_{\rm rel}$ for planetary events. Thus far, Penny et al. (2016) attempted
to compare the distance distribution of published microlensing planetary
systems with predictions based on a Galactic model and proposed a possibility
that the Galactic bulge might be devoid of planets. However, their results
were affected by an inhomogeneous sample and incorrect microlens parallax
measurements (Han et al., 2016).

In this letter, we compare the $\mu_{\rm rel}$ distribution for a given
$t_{\rm E}$ of 28 planetary events from the statistical sample by Suzuki et
al. (2016) with the distribution calculated using a Galactic model optimized
for use in microlensing studies (Koshimoto et al., 2021; Koshimoto & Ranc,
2021). We consider a power-law distribution, i.e., $P_{\rm host}\propto M_{\rm
L}^{m}R_{\rm L}^{r}$, for the planet-hosting probability for lens stars, where
$R_{\rm L}$ is the Galactocentric distance at which the lens is located. The
comparison of the data and model for various $(m,r)$ values enables us to
estimate for the first time the dependence of the planet-hosting probability
on the Galactocentric distance as $P_{\rm host}\propto R_{\rm L}^{0.2\pm 0.4}$
when a uniform prior for $m$ in $0<m<2$ is applied.

##  2 Method

This work aims to estimate the dependence of the planet-hosting probability on
the Galactic location by comparing the $\mu_{\rm rel}$ distribution observed
in planetary microlensing events with that predicted by a Galactic model.
Although the microlens parallax and/or lens brightness have already been
measured for some of the events in the sample, we here focus on the $t_{\rm
E}$ and $\mu_{\rm rel}$ distributions to avoid any bias caused by including
them. The $\mu_{\rm rel}$ distribution alone has no and little information
about the lens mass and distance, respectively (Penny et al., 2016). However,
when combined with $t_{\rm E}$, $\mu_{\rm rel}$ is equivalent to $\theta_{\rm
E}$, which yields a mass–distance relation and allows us to extract mass and
distance information.

Koshimoto & Bennett (2020) showed that an observed $\mu_{\rm rel}$
distribution can be compared with a model-expected distribution without
considering a detection efficiency correction once $t_{\rm E}$ is fixed, i.e.,

$\displaystyle f_{\rm obs}(\mu_{\rm rel}\,|\,t_{\rm E})\propto\Gamma_{\rm
Gal}(\mu_{\rm rel}\,|\,t_{\rm E}),$ (3)

where $f_{\rm obs}(\mu_{\rm rel}\,|\,t_{\rm E})$ and $\Gamma_{\rm
Gal}(\mu_{\rm rel}\,|\,t_{\rm E})$ are observed and expected $\mu_{\rm rel}$
distributions for given $t_{\rm E}$, respectively. $\Gamma_{\rm Gal}(\mu_{\rm
rel}\,|\,t_{\rm E})=\Gamma_{\rm Gal}(\mu_{\rm rel},t_{\rm E})/\Gamma_{\rm
Gal}(t_{\rm E})$, and $\Gamma_{\rm Gal}({\bm{x}})$ is the microlensing event
rate of events with parameter ${\bm{x}}$, which is calculated using the
Galactic model explained in Section 4. The observed distribution is given by
the sample described in Section 3.

A key idea behind Eq. (3) is that the detection efficiency of a microlensing
event depends on $t_{\rm E}$ but not on $\mu_{\rm rel}$; thus, detection
efficiency calculations are not necessary to compare observations with
simulations for a given $t_{\rm E}$. Technically, this is not true because the
detection of a planetary signal in a light curve depends on the source radius
crossing time, $t_{*}=\theta_{*}/\mu_{\rm rel}$, where $\theta_{*}$ is the
angular source radius. However, this is expected to have little effect on our
results because the dependence of the planet detection efficiency on $t_{*}$
is negligibly small for a mass-ratio of $q\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10^{-4}$ (Suzuki et al., 2016),
which dominates our sample.

##  3 Planetary Microlensing Event Sample

We used 28 planetary microlensing events, which combine 22 planetary events
detected by the MOA-II microlensing survey from 2007 to 2012 (Suzuki et al.,
2016) and 6 planetary events from Gould et al. (2010) and Cassan et al.
(2012). This is the same sample as the combined MOA+$\mu$FUN+PLANET sample in
Suzuki et al. (2016), except for OGLE-2011-BLG-0950 (Choi et al., 2012). The
event was excluded from our sample because both recent high-angular resolution
imaging observations by Terry et al. (in preparation) and the model-based
prior probability calculation by Koshimoto et al. (2021) suggest that
OGLE-2011-BLG-0950 was likely to be a stellar binary-lens event.

There are two events that have degenerate solutions with different $\mu_{\rm
rel}$ values in the sample. MOA-2011-BLG-262 (Bennett et al., 2014) has the
fast solution with $\mu_{\rm rel}=19.6\pm 1.6$ mas/yr and the slow solution
with $\mu_{\rm rel}=11.6\pm 0.9$ mas/yr. We use only the slow solution in our
analysis because it has a much larger prior probability as discussed in
Bennett et al. (2014). MOA-2010-BLG-328 (Furusawa et al., 2013) has two
parallax solutions ($u_{0}<0$ and $u_{0}>0$) and one xallarap solution, and
the $\mu_{\rm rel}$ values are measured as $5.71\pm 0.70~{}{\rm mas/yr}$,
$4.72\pm 0.79~{}{\rm mas/yr}$, and $4.03\pm 0.26~{}{\rm mas/yr}$,
respectively. In the likelihood calculation described in Section 5, we equally
combine the two parallax solutions, and then equally combine the combined
parallax solution with the xallarap solution, i.e., we use $\sqrt{{\cal
L}_{\rm xalla}\sqrt{{\cal L}_{\rm para+}\,{\cal L}_{\rm para-}}}$ as the
likelihood for this event.

The black open circles in Fig. 1 represent the $t_{\rm E}$ and $\mu_{\rm rel}$
values of the 28 planetary events; two of them, MOA-2007-BLG-192 (Bennett et
al., 2008) and MOA-2011-BLG-322 (Shvartzvald et al., 2014), have only lower
limit measurements on $\mu_{\rm rel}$. A more detailed description of Fig. 1
is given in Section 5.

Figure 1: Comparison of the $\mu_{\rm rel}$ values of the 28 planetary events
from the Suzuki et al. (2016) combined sample and Galactic model. The black
open circles show the data, while the color map shows the model-calculated
$\mu_{\rm rel}$ distribution for a given $t_{\rm E}$, $\Gamma_{\rm
Gal}(\mu_{\rm rel}\,|\,t_{\rm E})$, as a function of fixed $t_{\rm E}$. The
black solid, blue dashed, and magenta dotted lines indicate the median,
$1~{}\sigma$, and $2~{}\sigma$ for $\Gamma_{\rm Gal}(\mu_{\rm rel}\,|\,t_{\rm
E})$, respectively. (a) Model with $P_{\rm host}\propto{\rm const}$, i.e., all
stars are equally likely to host planets independent of their masses or
locations in our Galaxy. (b) Model with $P_{\rm host}\propto M_{\rm
L}^{0.2}R_{\rm L}^{0.2}$, corresponding to the best-fit grid that gives the
maximum likelihood.

##  4 Galactic Model for Planet-hosting Stars

To calculate the model of the $\mu_{\rm rel}$ distribution for a given $t_{\rm
E}$, $\Gamma_{\rm Gal}(\mu_{\rm rel}\,|\,t_{\rm E})$, we need a Galactic model
for planet-hosting stars. This combines a model of planet-hosting probability,
including the fit parameters investigated in this work, and a Galactic model
for random-selected stars, which refers to a combination of the stellar mass
function, stellar density distribution, and stellar velocity distribution in
our Galaxy.

We consider that the planet-hosting probability for the lens depends on the
lens mass $M_{\rm L}$ and Galactocentric distance $R_{\rm L}$ at which the
lens is located, where Galactocentric distance refers to the radius in the
cylindrical coordinate system with the Galactic center at the origin. The
exact relation remains unknown; thus, in this study, a power law was adopted,

$\displaystyle P_{\rm host}(M_{\rm L},R_{\rm L})\propto M_{\rm L}^{m}R_{\rm
L}^{r},$ (4)

where $m$ and $r$ are the fit parameters for the powers of $M_{\rm L}$ and
$R_{\rm L}$, respectively. Note that we used $P_{\rm host}(R_{\rm L})={\rm
const.}$ when $R_{\rm L}<50~{}{\rm pc}$ to avoid singularity with $r<0$. In
the Appendix, we consider a dichotomous model for $P_{\rm host}(R_{\rm L})$
that does not have this singularity.

###  4.1 Galactic Model for Random-selected Stars

Koshimoto et al. (2021) developed a parametric Galactic model by fitting to
Gaia DR2 velocity data (Gaia Collaboration et al., 2018), OGLE-III red clump
star count data (Nataf et al., 2013), VIRAC proper motion data (Smith et al.,
2018; Clarke et al., 2019), BRAVA radial velocity data (Rich et al., 2007;
Kunder et al., 2012), and OGLE-IV star count and microlensing rate data (Mróz
et al., 2017, 2019). All the data except for the Gaia DR2 data correspond to a
bulge region of the sky, and the models are optimized for use in microlensing
studies toward the Galactic bulge.

There are four versions of models developed in Koshimoto et al. (2021), which
are denoted by E, G, E+EX, and G+GX. Each model consists of a multicomponent
thin disk, thick disk, and barred bulge. For the bulge density profile,
exponential functions are used in the E and E+EX models, whereas Gaussian
functions are used in the G and G+GX models. The E and G models have a single-
component bulge density profile, whereas the E+EX and G+GX models have an
additional component to represent the X-shape structure in the bulge
(McWilliam & Zoccali, 2010; Nataf et al., 2010). The different bulge profiles
used in the four models led to different best-fit parameters for the stellar
density, velocity, and initial mass function for all four models.

In this study, we used the E+EX model as our fiducial Galactic model for
random-selected stars because it was the most consistent with the data used in
Koshimoto et al. (2021). Nevertheless, we also used the other three models to
evaluate the systematic errors in our measurement of the planet-hosting
probability dependence.

###  4.2 Possibilities of Planets around White Dwarfs or Close-binaries

Our sample includes at least one circumbinary planet (OGLE-2007-BLG-349,
Bennett et al., 2016). Because the number of events with lens brightness
measurements in our sample is limited, there is a non-negligible chance that
lens objects of other events may also consist of circumbinary planets or even
planets around white dwarf hosts. Thus, both possibilities should be
considered in the calculation of the event rate $\Gamma_{\rm Gal}(\mu_{\rm
rel}\,|\,t_{\rm E})$.

We introduce two parameters to account for planet frequencies relative to
single stars: $f_{\rm WD}$ for planets around white dwarfs and $f_{\rm CB}$
for planets around tight close binaries. Values of $f_{\rm WD}=1$ or $f_{\rm
CB}=1$ indicate the same planet frequency as for a single star with the same
mass; if $f_{\rm WD}=0$ or $f_{\rm CB}=0$, there is zero probability of
hosting planets. Determining the details of the planet abundance around white
dwarfs or close binary systems is beyond the scope of this study; thus, we
adopted ($f_{\rm WD}$, $f_{\rm CB}$) = (1, 1) for our fiducial model and used
cases ($f_{\rm WD}$, $f_{\rm CB}$) = (0, 1), (1, 0), and (0, 0) to evaluate
the systematic errors in our $r$ estimate.

To consider the white dwarf population and tight close-binary population in
the model, we followed Koshimoto et al. (2021), where the initial-final mass
relation of $M_{\rm WD}=0.109M_{\rm ini}+0.394M_{\odot}$ (Kalirai et al.,
2008) is used for the white dwarfs. Further, the binary distribution developed
by Koshimoto et al. (2020), together with a detection threshold based on a
combination of the central caustic size of a hypothetical binary and the event
impact parameter $u_{0}$, are used for tight close-binaries.

We assume that the planet-hosting probability for a neutron star or black hole
is zero. The assumption is based on the study by Behrens et al. (2020) that
found no planets in an 11-year dataset for 45 pulsars. If this rareness of
planets around neutron stars is attributed to the progenitor’s explosion
(i.e., supernova), planets are also likely to be rare around black holes.

##  5 Maximum Likelihood Analysis

The color maps in Fig. 1 are two examples of the model-calculated $\mu_{\rm
rel}$ distribution for given $t_{\rm E}$, $\Gamma_{\rm Gal}(\mu_{\rm
rel}\,|\,t_{\rm E})$, $(m,r)=(0,0)$ in Fig. 1(a) and $(m,r)=(0.2,0.2)$ in Fig.
1(b). The distributions are calculated over the range $0.50<\log(t_{\rm
E}/{\rm days})<2.20$. First, this range is divided into 34 bins of width 0.05
dex each. Then, $10^{5}$ artificial events are generated using our fiducial
Galactic model for planet-hosting stars, i.e., the E+EX model with ($f_{\rm
WD}$, $f_{\rm CB}$) = (1, 1). We selected a Galactic coordinate of
$(l,b)=(1.0^{\circ},-2.2^{\circ})$ for these plots. The coordinates of each
event are used in the relative likelihood calculations presented below.

Fig. 1 (a) shows that the observed $\mu_{\rm rel}$ distribution is already in
good agreement with the model, in which no planet-hosting probability
dependence on the stellar mass or position in our Galaxy is assumed. This
means that $(m,r)=(0,0)$ is an acceptable planet-hosting probability
dependence under the current constraint imposed by the data. Below, we
calculate the relative likelihood for various combinations of $(m,r)$ and
obtain the probability distributions for these parameters.

###  5.1 Definition of Likelihood

We define the likelihood for a combination of $m$ and $r$ as

$\displaystyle{\cal L}(m,r)=\prod_{i}\Gamma_{\rm Gal}(\mu_{{\rm
rel},i}\,|\,t_{{\rm E},i}\,;\,m,r),$ (5)

where $\mu_{{\rm rel},i}$ and $t_{{\rm E},i}$ are the ones observed for $i$th
event, and $\Gamma_{\rm Gal}(\mu_{{\rm rel},i}\,|\,t_{{\rm E},i}\,;\,m,r)$ is
the model-calculated probability of $\mu_{{\rm rel},i}$ for given $t_{{\rm
E},i}$ when the planet-hosting probability is $P_{\rm host}\propto M_{\rm
L}^{m}R_{\rm L}^{r}$. The product is taken over all 28 events in our sample.

For each event $i$, in the calculation of $\Gamma_{\rm Gal}(\mu_{{\rm
rel},i}\,|\,t_{{\rm E},i}\,;\,m,r)$, the corresponding parameters are used,
e.g., the Galactic coordinate and impact parameter $u_{0}$ affecting the
close-binary possibilities.

###  5.2 Results

Figure 2: Relative likelihood distribution as a function of $m$ and $r$ for
the E+EX fiducial Galactic model with ($f_{\rm WD}$, $f_{\rm CB}$) = (1, 1).
In the top panel, probability distributions of $r$ are shown in solid black
and dotted red, integrated ${\cal L}(m,r)$ over $-2<m<2$ and $0<m<2$,
respectively. The dark and light gray lines in the edge panels are the
marginalized probability distributions when ($f_{\rm WD}$, $f_{\rm CB}$) = (0,
1) and ($f_{\rm WD}$, $f_{\rm CB}$) = (1, 0), respectively. (Note that the
light gray line almost overlaps with the black line.)

We calculated the likelihood given by Eq. (5) over the ranges of $-2<m<2$ and
$-2<r<2$ at intervals of 0.2, for both $m$ and $r$. This corresponds to
applying uniform prior distributions for the two parameters between -2 and 2,
which reflects that there are moderate fractions of planets in all categories
of low-mass lenses (e.g., MOA-2008-BLG-310, Janczak et al., 2010; Bhattacharya
et al., 2017; Koshimoto et al., 2020), high-mass lens (e.g.,
OGLE-2012-BLG-0026, Han et al., 2013; Beaulieu et al., 2016), distant lens
(e.g., MOA-2007-BLG-400, Dong et al., 2009; Bhattacharya et al., 2021), and
close lens (e.g., OGLE-2006-BLG-109, Gaudi et al., 2008; Bennett et al.,
2010). The planet frequency dependence on $M_{\rm L}$ or $R_{\rm L}$ would not
be as extreme as $|m|>2$ or $|r|>2$.

Fig. 2 shows the results, where a clear negative correlation between $m$ and
$r$ is seen. The shape of the correlation is attributed to the mass–distance
relation given by the angular Einstein radius $\theta_{\rm E}$, a parameter
equivalent to $\mu_{\rm rel}$ for a given $t_{\rm E}$. Because certain
$\theta_{\rm E}$ values can be explained by both a low-mass lens close to the
Sun and a high-mass lens close to the source, the data distribution can be
reproduced by both small $m$ with large $r$ and large $m$ with small $r$. This
degeneracy can be disentangled by including observational information of
either the microlens parallax, lens brightness, or Einstein radius crossing
time $t_{\rm E}$ with a proper detection efficiency correction, which all are
interesting future works.

On the other hand, the location of the negatively correlated distribution on
the $r$ vs. $m$ plane is determined by the data and Galactic model. Fig. 2
shows a high likelihood value at $(m,r)=(0,0)$, as expected from Fig. 1 (a).
For comparison, we show the $\Gamma_{\rm Gal}(\mu_{\rm rel}\,|\,t_{\rm E})$
distribution using the best-fit grid value of $(m,r)=(0.2,0.2)$ in Fig. 1 (b).
The correlation between $r$ and $m$ becomes weaker when $m$ increases, and a
stronger constraint can be placed on $r$ than $m$, for which we have estimated
$r=0.4^{+0.6}_{-0.4}$.

Because it is difficult to put a useful constraint on the $m$ value based on
the current data, we considered a stronger uniform prior on $m$ in the range
$0<m<2$. This is a plausible prior considering the ongoing high-angular
resolution follow-up program on past events (Bennett, 2018), which revealed
that two events (out of seven published) have lens host stars that are much
more massive than those expected by Bayesian analysis using a Galactic model
for random-selected stars (Vandorou et al., 2020; Bhattacharya et al.,
2021).222 Our preliminary result from a statistical study on the Suzuki et al.
(2016) sample based on mass measurements implies $m>0$ with a $3~{}\sigma$
confidence level, as reported by DPB in the Exoplanet Demographics Conference
in November 2020. The talk can be found on
https://www.youtube.com/watch?v=5TpOKjHSS10 (catalog
https://www.youtube.com/watch?v=5TpOKjHSS10). For giant planets, the
possibility of $m>0$ is also supported by other techniques. Johnson et al.
(2010) estimated $m=1.0\pm 0.3$ in $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2.5~{}$AU from a radial velocity
survey, and Nielsen et al. (2019) estimated $m=2.0\pm 1.0$ in 10–100 AU from a
direct imaging survey.

As shown by the red dotted curve in the top panel of Fig. 2, applying the
narrower prior of $0<m<2$ can almost entirely exclude the possibility of
$r>1$. The median and $1~{}\sigma$ error is $r=0.2\pm 0.4$ for the fiducial
E+EX Galactic model with ($f_{\rm WD}$, $f_{\rm CB}$) = (1, 1).

Finally, to evaluate the systematic error due to model selection, we repeated
this analysis using the E, G, and G+GX models and for ($f_{\rm WD}$, $f_{\rm
CB}$) = (1, 0), (0, 1), and (0, 0). We found that the median value of $r$
becomes 0.3–0.5 and 0.2–0.3 for the uniform priors in $-2<m<2$ and $0<m<2$,
respectively. Thus, the systematic errors seem to be much smaller than the
statistical errors. Table 1 summarizes the results of our likelihood analysis.

Table 1: Results of the maximum likelihood analysis.  Range of uniform prior | $r$ |  |  $m$ ccListed for completeness. The estimates of $m$ are dominated by the prior applied.
---|---|---|---
| Fiducial aaMedian and $1~{}\sigma$ error for the E+EX model and $(f_{\rm WD},f_{\rm CB})=(1,1)$. | Sys. range bbVariation of median values when other Galactic models or other combinations of $(f_{\rm WD},f_{\rm CB})$ are applied. |  | Fiducial aaMedian and $1~{}\sigma$ error for the E+EX model and $(f_{\rm WD},f_{\rm CB})=(1,1)$.
$-2<m<2$ | $0.4^{+0.6}_{-0.4}$ | 0.3–0.5 |  | $0.2\pm 1.0$
$0<m<2$ | $0.2\pm 0.4$ | 0.2–0.3 |  | $0.7^{+0.8}_{-0.6}$

##  6 Discussion and Conclusion

We estimated the planet-hosting probability dependence of $P_{\rm host}\propto
R_{\rm L}^{0.2\pm 0.4}$ under a uniform prior distribution of $m$ in $0<m<2$,
which suggests no large dependence of the planet frequency on the
Galactocentric distance. The bulge region has a very different stellar
environment from the solar neighborhood, including stellar densities $\lower
2.15277pt\hbox{$\;\buildrel>\over{\sim}\;$}10$ times higher and an older,
alpha-enhanced stellar population. Observations of the solar neighborhood have
shown that there is a correlation between stellar metallicity and the
occurrence of giant planets (Fischer & Valenti, 2005; Johnson et al., 2010),
and have also suggested that close encounters with other stars may affect the
evolution of planetary systems (Winter et al., 2020). Therefore, due to the
abovementioned environmental differences, the planet frequency in the bulge
may differ significantly from that in the solar neighborhood. Although our
results are still inconclusive, they might imply that cold planets orbiting
beyond the H2O snow line also commonly exist in the bulge region regardless of
such differences.

Because all exoplanet detection techniques have sensitivities that complement
each other, it is important to combine different techniques for a
comprehensive understanding of planet formation and evolution processes (Gaudi
et al., 2020). One of the largest uncertainties when comparing the exoplanet
population discovered via other exoplanet detection methods with microlensing
planets is a possible large dependence of the planet frequency on Galactic
location. Our results show that such a dependence is not very large, and one
might be able to compare them without considering the difference in Galactic
location. A small dependence of the planet frequency on the Galactic location
is also supported by Suzuki et al. (2016). Suzuki et al. (2016) showed that
the planet frequency from the MOA-II microlensing survey—which is an averaged
value for stars in the galactic disk and bulge—is consistent with the
frequency of cold gas giants from radial velocity studies (Bonfils et al.,
2013; Montet et al., 2014).

On the other hand, the full consistency with $r=0$ contrasts with the analysis
by Penny et al. (2016), who found a small $p$-value of 5.0 $\times 10^{-4}$
for a model in which the planetary frequency of the bulge is the same as that
of the disk. The difference is most likely due to the contamination of
excessively close lenses with incorrect parallax measurements in their sample.
Penny et al. (2016) themselves discussed this possibility, and questioned the
results of several events with large microlens parallax values. In fact, among
the questioned events, OGLE-2013-BLG-0723 was a stellar binary event (Han et
al., 2016). The large microlens parallax claimed for MOA-2007-BLG-192 (Bennett
et al., 2008) has recently been shown to disappear when the data are detrended
for color-dependent differential refraction (Bennett et al., 2012), although
this has not yet been published. Nevertheless, our results are still
consistent with the Penny et al. (2016)’s estimate of the bulge-to-disk ratio
of planet frequency, $f_{\rm bulge}<0.54$, as shown in the Appendix.

Our estimate of $r=0.2\pm 0.4$ still has a moderately large uncertainty, and
there could be a non-negligible dependence of planet frequency on the
Galactocentric distance $R_{\rm L}$. In particular, although $m>0$ seems to be
thus far plausible (Vandorou et al., 2020; Bhattacharya et al., 2021), if the
possibility of $m<0$ is considered, the uncertainty and preference for disk
planets could increase. This is owing to the negative correlation between $m$
and $r$, as shown in Fig. 2, which is attributed to our use of only one
parameter (i.e., $\theta_{\rm E}$) that provides a lens mass–distance relation
in the analysis. Thus, further constraints should emerge by including
constraints from other mass–distance relations, which can be provided by a
statistical sample of either microlens parallax or lens brightness
measurements.

NK was supported by the JSPS overseas research fellowship. DPB and NK were
supported by NASA through grant NASA-80NSSC18K0274 and award number
80GSFC17M0002. DS was supported by JSPS KAKENHI Grant Number JP19KK082 and
JP20H04754.

##  Appendix A Dichotomous Model for Planet-hosting Probability

One might think that the power law model, $P_{\rm host}\propto R_{\rm L}^{r}$,
is not physical because it diverges at $R_{\rm L}=0$ when $r<0$. Here, we
consider a dichotomous model for the dependence of the planet-hosting
probability on the lens Galactic location,

$P_{\rm host}(R_{\rm L})=\begin{cases}P_{\rm host,D}&\text{ when $R_{\rm
L}>2{\rm kpc}$}\\\ P_{\rm host,B}&\text{ when $R_{\rm L}\leq 2{\rm
kpc}$,}\end{cases}$ (A1)

and we use the ratio between the two constants, $f_{\rm B/D}\equiv P_{\rm
host,B}/P_{\rm host,D}$, as the fit parameter instead of $r$. This is a more
directly comparable model to that used by Penny et al. (2016), who considered
the bulge-to-disk ratio of planet frequency in the Han & Gould (2003) Galactic
model, $f_{\rm bulge}$. Note that there is still a subtle difference between
the two parameters—$f_{\rm B/D}=0$ indicates that there are no planets at
$R_{\rm L}\leq 2{\rm kpc}$, whereas $f_{\rm bulge}=0$ allows planet-hosting
stars in the disk component to exist at $R_{\rm L}\leq 2{\rm kpc}$.

Figure 3: Same as Fig. 2, but for $f_{\rm B/D}$ instead of $r$. Table 2: $f_{\rm B/D}$ estimated from the maximum likelihood analysis. Range of uniform prior | $f_{\rm B/D}$
---|---
| Fiducial aaMedian and $1~{}\sigma$ error with the E+EX model and $(f_{\rm WD},f_{\rm CB})=(1,1)$. | Sys. range bbVariation of median values when other Galactic models or other combinations of $(f_{\rm WD},f_{\rm CB})$ are applied.
$-2<m<2$ | $0.48^{+0.59}_{-0.33}$ | 0.42 – 0.52
$0<m<2$ | $0.63^{+0.61}_{-0.33}$ | 0.56 – 0.65

Fig. 3 and Table 2 show the results of the maximum likelihood analysis using
$f_{\rm B/D}$ instead of $r$, where the uniform prior for $\log f_{\rm B/D}$
in $-2<\log f_{\rm B/D}<2$ is assumed. With the uniform prior for $m$ in
$0<m<2$, we estimate $f_{\rm B/D}=0.63^{+0.61}_{-0.33}$, which is consistent
with both $f_{\rm B/D}=1$ and $f_{\rm bulge}<0.54$; the upper limit estimate
for $f_{\rm bulge}$ was taken from a $p$-value threshold of 0.01 by Penny et
al. (2016).

## References

  * Akeson et al. (2013) Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989. doi:10.1086/672273
  * Beaulieu et al. (2016) Beaulieu, J.-P., Bennett, D. P., Batista, V., et al. 2016, ApJ, 824, 83. doi:10.3847/0004-637X/824/2/83
  * Behrens et al. (2020) Behrens, E. A., Ransom, S. M., Madison, D. R., et al. 2020, ApJ, 893, L8. doi:10.3847/2041-8213/ab8121
  * Bennett (2018) Bennett, D. 2018, Keck Observatory Archive N139, 172
  * Bennett et al. (2014) Bennett, D. P., Batista, V., Bond, I. A., et al. 2014, ApJ, 785, 155. doi:10.1088/0004-637X/785/2/155
  * Bennett et al. (2008) Bennett, D. P., Bond, I. A., Udalski, A., et al. 2008, ApJ, 684, 663. doi:10.1086/589940
  * Bennett et al. (2010) Bennett, D. P., Rhie, S. H., Nikolaev, S., et al. 2010, ApJ, 713, 837. doi:10.1088/0004-637X/713/2/837
  * Bennett et al. (2016) Bennett, D. P., Rhie, S. H., Udalski, A., et al. 2016, AJ, 152, 125. doi:10.3847/0004-6256/152/5/125
  * Bennett et al. (2012) Bennett, D. P., Sumi, T., Bond, I. A., et al. 2012, ApJ, 757, 119. doi:10.1088/0004-637X/757/2/119
  * Bhattacharya et al. (2018) Bhattacharya, A., Beaulieu, J.-P., Bennett, D. P., et al. 2018, AJ, 156, 289
  * Bhattacharya et al. (2017) Bhattacharya, A., Bennett, D. P., Anderson, J., et al. 2017, AJ, 154, 59. doi:10.3847/1538-3881/aa7b80
  * Bhattacharya et al. (2021) Bhattacharya, A., Bennett, D. P., Beaulieu, J. P., et al. 2021, arXiv:2009.02329, accepted
  * Bonfils et al. (2013) Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109. doi:10.1051/0004-6361/201014704
  * Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167. doi:10.1038/nature10684
  * Choi et al. (2012) Choi, J.-Y., Shin, I.-G., Han, C., et al. 2012, ApJ, 756, 48. doi:10.1088/0004-637X/756/1/48
  * Clarke et al. (2019) Clarke, J. P., Wegg, C., Gerhard, O., et al. 2019, MNRAS, 489, 3519. doi:10.1093/mnras/stz2382
  * Dong et al. (2009) Dong, S., Bond, I. A., Gould, A., et al. 2009, ApJ, 698, 1826. doi:10.1088/0004-637X/698/2/1826
  * Fischer & Valenti (2005) Fischer, D. A. & Valenti, J. 2005, ApJ, 622, 1102. doi:10.1086/428383
  * Furusawa et al. (2013) Furusawa, K., Udalski, A., Sumi, T., et al. 2013, ApJ, 779, 91. doi:10.1088/0004-637X/779/2/91
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 616, A11
  * Gaudi (2012) Gaudi, B. S. 2012, ARA&A, 50, 411. doi:10.1146/annurev-astro-081811-125518
  * Gaudi et al. (2008) Gaudi, B. S., Bennett, D. P., Udalski, A., et al. 2008, Science, 319, 927. doi:10.1126/science.1151947
  * Gaudi et al. (2020) Gaudi, B. S., Christiansen, J. L., & Meyer, M. R. 2020, arXiv:2011.04703
  * Gould et al. (2010) Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073. doi:10.1088/0004-637X/720/2/1073
  * Han et al. (2016) Han, C., Bennett, D. P., Udalski, A., et al. 2016, ApJ, 825, 8. doi:10.3847/0004-637X/825/1/8
  * Han & Gould (2003) Han, C. & Gould, A. 2003, ApJ, 592, 172. doi:10.1086/375706
  * Han et al. (2013) Han, C., Udalski, A., Choi, J.-Y., et al. 2013, ApJ, 762, L28. doi:10.1088/2041-8205/762/2/L28
  * Janczak et al. (2010) Janczak, J., Fukui, A., Dong, S., et al. 2010, ApJ, 711, 731. doi:10.1088/0004-637X/711/2/731
  * Johnson et al. (2010) Johnson, J. A., Aller, K. M., Howard, A. W., et al. 2010, PASP, 122, 905. doi:10.1086/655775
  * Kalirai et al. (2008) Kalirai, J. S., Hansen, B. M. S., Kelson, D. D., et al. 2008, ApJ, 676, 594. doi:10.1086/527028
  * Koshimoto et al. (2021) Koshimoto, N., Baba, J., & Bennett, D. P. 2021, ApJ, 917, 78. doi:10.3847/1538-4357/ac07a8
  * Koshimoto & Ranc (2021) Koshimoto, N. & Ranc, C. 2021, Zenodo.4784948. doi: 10.5281/zenodo.4784948
  * Koshimoto & Bennett (2020) Koshimoto, N. & Bennett, D. P. 2020, AJ, 160, 177. doi:10.3847/1538-3881/abaf4e
  * Koshimoto et al. (2020) Koshimoto, N., Bennett, D. P., & Suzuki, D. 2020, AJ, 159, 268. doi:10.3847/1538-3881/ab8adf
  * Kunder et al. (2012) Kunder, A., Koch, A., Rich, R. M., et al. 2012, AJ, 143, 57. doi:10.1088/0004-6256/143/3/57
  * McWilliam & Zoccali (2010) McWilliam, A. & Zoccali, M. 2010, ApJ, 724, 1491. doi:10.1088/0004-637X/724/2/1491
  * Montet et al. (2014) Montet, B. T., Crepp, J. R., Johnson, J. A., et al. 2014, ApJ, 781, 28. doi:10.1088/0004-637X/781/1/28
  * Mróz et al. (2017) Mróz, P., Udalski, A., Skowron, J., et al. 2017, Nature, 548, 183
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, J., et al. 2019, ApJS, 244, 29. doi:10.3847/1538-4365/ab426b
  * Nataf et al. (2013) Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88. doi:10.1088/0004-637X/769/2/88
  * Nataf et al. (2010) Nataf, D. M., Udalski, A., Gould, A., et al. 2010, ApJ, 721, L28. doi:10.1088/2041-8205/721/1/L28
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13. doi:10.3847/1538-3881/ab16e9
  * Penny et al. (2016) Penny, M. T., Henderson, C. B., & Clanton, C. 2016, ApJ, 830, 150
  * Rich et al. (2007) Rich, R. M., Reitzel, D. B., Howard, C. D., et al. 2007, ApJ, 658, L29. doi:10.1086/513509
  * Shvartzvald et al. (2014) Shvartzvald, Y., Maoz, D., Kaspi, S., et al. 2014, MNRAS, 439, 604. doi:10.1093/mnras/stt2477
  * Smith et al. (2018) Smith, L. C., Lucas, P. W., Kurtev, R., et al. 2018, MNRAS, 474, 1826. doi:10.1093/mnras/stx2789
  * Suzuki et al. (2016) Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145
  * Vandorou et al. (2020) Vandorou, A., Bennett, D. P., Beaulieu, J.-P., et al. 2020, AJ, 160, 121. doi:10.3847/1538-3881/aba2d3
  * Winter et al. (2020) Winter, A. J., Kruijssen, J. M. D., Longmore, S. N., et al. 2020, Nature, 586, 528. doi:10.1038/s41586-020-2800-0


# Counting cusp excursions of reciprocal geodesics

Ara Basmajian  PhD Program in Mathematics, The Graduate Center, CUNY, 365
Fifth ave., NY,NY, 10016 and Department of Mathematics and Statistics, Hunter
College, CUNY, 695 Park ave. NY, NY, 10065, USA<EMAIL_ADDRESS>and
Robert Suzzi Valli  Department of Mathematics, Manhattan College, 4513
Manhattan College Pkwy., Riverdale, NY, 10471<EMAIL_ADDRESS>Dedicated to Ravi S. Kulkarni on his 80th birthday

###### Abstract.

For a fixed cusp neighborhood (determined by depth $D$) of the modular
surface, we investigate the class of reciprocal geodesics that enter this
neighborhood (called a cusp excursion) a fixed number of times. Since
reciprocal geodesics traverse their image an even number of times such a
geodesic only has an even number of cusp excursions, say $2n$. Denoting this
class of geodesics as $\mathcal{E}_{2n}^{D}$, we derive several growth rate
results for $|\\{\gamma\in\mathcal{E}_{2n}^{D}:|\gamma|=4t\\}|$, where
$|\gamma|$ is the word length of $\gamma$ in
$\mathbb{Z}_{2}\ast\mathbb{Z}_{3}$ with respect to the generators of the
factors.

###### Key words and phrases:

Reciprocal geodesic, modular group, modular surface, cusp excursion, word
length

###### 1991 Mathematics Subject Classification:

Primary 20F69, 32G15, 57K20; Secondary 20H10, 53C22.

The first author was supported by a PSC-CUNY Grant and a grant from the Simons
foundation (359956, A.B.)

##  1\. Introduction

For a fixed cusp neighborhood determined by its depth $D$, the paper [2] was
interested in counting the number of so called $D$-low lying reciprocal
geodesics on the modular surface. That is, in counting the reciprocal
geodesics that stay in the compact subsurface determined by chopping off the
$D$-cusp neighborhood. In this paper, we are interested in counting the number
of reciprocal geodesics on the modular surface that leave this compact
subsurface, called a cusp excursion, a fixed number of times. Since reciprocal
geodesics traverse their image an even number of times such a geodesic only
has an even number of cusp excursions, say $2n$. Denoting this class of
geodesics as $\mathcal{E}_{2n}^{D}$, we derive growth rate results (Theorem
1.1) for $|\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}|$, and
$|\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|=4t\\}|$. Before we state our
theorem we first define our terms and set up language.

Consider the group $G=\mathbb{Z}_{2}\ast\mathbb{Z}_{3}$ with its natural
symmetric generators $a$, $b$, and $b^{-1}$ where $a\in\mathbb{Z}_{2}$ has
order two, and $b\in\mathbb{Z}_{3}$ has order 3. The representation of $G$
given by, $a\mapsto A$ and $b\mapsto B$ where,

$A=\begin{pmatrix}0&-1\\\ 1&0\end{pmatrix},\text{ and
}B=\begin{pmatrix}1&-1\\\ 1&0\end{pmatrix}$

is a discrete faithful representation with image the modular group
$PSL(2,\mathbb{Z}).$ The quotient $S=\mathbb{H}/PSL(2,\mathbb{Z})$, called the
modular surface, is a generalized pair of pants of genus zero and signature
$(2,3,\infty)$. For $\gamma$ a closed geodesic on $S$, its word length,
denoted $|\gamma|$, is the minimal word length (with respect to the natural
generators of the factors in $\mathbb{Z}_{2}\ast\mathbb{Z}_{3}$) among the
elements of the conjugacy class of a lift in $PSL(2,\mathbb{Z})$.

Now, it is well-known that if an orbifold surface has a cusp then it has a
cusp neighborhood, $\text{say }\mathcal{C}$, of area one with boundary a
horocycle segment of length one. For $D$ a positive integer, we define the
$D$-thick part of $S$, denoted $S_{D}$, to be the smallest compact subsurface
of $S$ with horocycle boundary which contains all the closed geodesics which
when written in normal form, $ab^{\epsilon_{1}}...ab^{\epsilon_{t}}$,
($\epsilon_{i}=\pm 1,\epsilon_{1}\neq\epsilon_{t}$) have at most $D$
consecutive $\epsilon_{i}^{\prime}s$ of the same sign. The number of
consecutive $\epsilon_{i}^{\prime}s$ of the same sign (called a run)
correspond to how deep the geodesic goes into the cusp. (See for example Lemma
7.1 in [2]). Hence it is natural to associate the geometric depth into the
cusp with the size of a run; that is, the highest power of the parabolic. The
$D$-thick parts form a compact exhaustion of $S$.

We say that a closed geodesic (including reciprocal geodesics) $\gamma$ has
$k$ (cusp) excursions of depth greater than $D$ if during one periodic orbit
$\gamma$ enters $S-S_{D}$ exactly $k$-times. See Figure 1 for an example of a
cusp excursion where $k=1$.

A reciprocal geodesic on the modular surface is a geodesic that begins and
ends at the order two cone point, traversing its image twice if it is
primitive, and more generally traversing the image $2m$-times if it is the
$m^{th}$ power of a reciprocal geodesic. Its lift is the conjugacy class in
$PSL(2,\mathbb{Z})$ of a hyperbolic element with axis passing through an order
two fixed point. In the paper [2], reciprocal geodesics that stay in $S_{D}$
(so called $D$-low lying reciprocal geodesics) were studied. In this paper, we
are interested in the reciprocal geodesics that leave the $D$-thick part a
prescribed (necessarily even) number of times (cusp excursions). To that end,
for each integer $n\geq 0$, define

$\mathcal{E}_{2n}^{D}=\\{\gamma:\gamma\text{ a reciprocal geodesic that enters
$S-S_{D}$ exactly $2n$-times}\\}$

Equivalently, $\gamma\in\mathcal{E}_{2n}^{D}$ if when represented in normal
form its associated run sequence has exactly $2n$ parts each of size at least
$D+1$. With this in mind, the reciprocal geodesics are partitioned into the
disjoint sets

${{\bigcup}_{n=0}^{\infty}}\mathcal{E}_{2n}^{D}$

The set $\mathcal{E}_{0}^{D}$ is exactly the $D$-low lying reciprocals whose
asymptotic growth rate is $d_{D}\alpha_{D}^{t}$, where $d_{D}$ is given by the
coefficient in expression (2.1), and $\alpha_{D}$ is the unique positive root
of $z^{D}-z^{D-1}-\dots-1$. In stark contrast, the set of all reciprocal
geodesics has asymptotic growth $2^{t-1}$ (see [2] for details).

\begin{overpic}[scale={.45}]{CuspExcursion} \put(0.0,36.0){$3$}
\put(9.5,11.0){$2$} \put(42.5,6.5){$S_{D}$} \put(33.0,20.0){$\gamma$}
\end{overpic} Figure 1. A cusp excursion of depth greater than $D$ on $S$

So the first non-trivial excursions occur in $\mathcal{E}_{2}^{D}$. In this
note we study the growth of reciprocal geodesics with a fixed number of
excursions. Specifically, we are interested in how the word lengths of the
geodesics in $\\{\gamma\in\mathcal{E}_{2n}^{D}:|\gamma|=4t\\}$ grow as a
function of $t$. In this paper, we address the growth of reciprocal geodesics
with any number of excursions greater than depth $1$, and the growth of
reciprocal geodesics with two excursions of depth greater than $D$.

######  Theorem 1.1.

For any integers $n\geq 0$,

(1.1)
$\frac{\big{|}\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|=4t\\}\big{|}}{t^{2n}}\xrightarrow{t\to\infty}\frac{1}{2n!}$

For any integer $D\geq 2$,

(1.2)
$\frac{\big{|}\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}\big{|}}{t\alpha_{D}^{t}}\xrightarrow{t\to\infty}\frac{d_{D}^{2}}{\alpha_{D}^{D}(\alpha_{D}-1)}$

where $d_{D}=\frac{\alpha_{D}-1}{2+(D+1)(\alpha_{D}-2)}$, and $\alpha_{D}$ is
the unique positive root of $z^{D}-z^{D-1}-\dots-1$.

In the sequel, Theorem 1.1 and its proof are done in two parts, Theorems 3.1
and 3.3. For the set $\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|=4t\\}$, we in
fact obtain a precise formula, see Theorem 3.1. Moreover in the course of the
proof of Theorem 3.3, we derive coarse bounds (in closed form) for the set
$\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}.$ As a result, in both cases,
it is straightforward to write down the growth rates of
$|\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|\leq 4t\\}|$ and
$|\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|\leq 4t\\}|$, as
$t\rightarrow\infty$. Finally, we note that the asymptotic constant in
expression (1.2) is only dependent on the depth $D$ and goes to zero as
$D\rightarrow\infty$. The more general setting of counting any number of
excursions of depth bigger than $D$ will be addressed elsewhere.

The first author would like to thank the University of Michigan and the
University of Luxembourg for their support while parts of this paper were
being completed. In particular, thanks go to Richard D. Canary (University of
Michigan) and Hugo Parlier (University of Luxembourg) for their hospitality
during the respective visits.

##  2\. Some background and preliminaries

There is an extensive history surrounding the modular group as well as an
extensive history around the question of geodesic growth in a hyperbolic
surface. For the study of growth from a geometric point of view we refer to
the papers, [3, 4, 5, 6, 10, 12, 27, 33]. From an algebraic point of view for
the study of surface groups and free groups we refer to the papers, [7, 9, 19,
25, 29, 32, 37]. Papers involving normal forms, enumeration schemes for
curves, Farey arithmetic, and generating elements in a non-abelian free group
can be found in [1, 13, 14, 15, 16, 17, 18, 34]. Papers specifically about
reciprocal geodesics appear in [2, 3, 5, 11, 25, 33]. There is an extensive
literature on cusp excursions by a random geodesic on a hyperbolic surface
[20, 21, 22, 23, 24, 26, 28, 30, 31, 35, 36, 38]. In particular, the papers
[36, 20, 21, 22, 23, 24] investigate the relation between depth, return time,
and other invariants in various contexts including connections to number
theory.

A closed geodesic in $S$ corresponds to a conjugacy class of hyperbolic
elements in $G$. Using the fact that $ab$ and $ab^{-1}$ are parabolic
elements, a conjugacy class of hyperbolic elements in $PSL(2,\mathbb{Z})$ has
a representative of the form $ab^{\epsilon_{1}}...ab^{\epsilon_{t}}$, where
$\epsilon_{i}=\pm 1$ and $\epsilon_{1}\neq\epsilon_{t}$. Being a product of
the parabolic elements $ab$ and $ab^{-1}$, if $k$ consecutive $\epsilon_{i}$
have the same sign, then $\gamma$ wraps around the cusp $k$-times, and hence
the hyperbolic geometry dictates that the geodesic $\gamma$ goes deeper into
the cusp.

Now, to a $t$-tuple $(\epsilon_{1},...,\epsilon_{t})$ we associate a run
sequence (composition of $t$) $(r_{1},...,r_{k})$, where $r_{1}$ is the number
of consecutive $\epsilon_{i}$ of the same sign starting with $\epsilon_{1}$,
$r_{2}$ is the length of the next such sequence, and so on. For example,
$(1,1,1,-1,1,-1,-1,1)$ corresponds to the composition of 8, $(3,1,1,2,1)$.
Thus a composition of $t$, denoted $C_{t}$, is a sequence of positive integers
that sum to $t$. The individual summands of a composition are called its
parts. When the parts are at most $m$, we denote that subset of $C_{t}$ by
$C_{t,m}$. Lastly, we denote a composition with exactly $n$ parts that are
bigger than $D$, and all other parts at most $D$ by $C_{t}^{n,D}$.

Geometrically the length of the run corresponds to how deep the geodesic goes
into the cusp. (See lemma 7.1 in [2]). Namely, for $D>0$, we say a closed
geodesic (reciprocal geodesics are included among the closed geodesics) is
said to have a cusp excursion (excursion for short) of depth bigger than $D$
if there is a run longer than $D$ in its normal form. For reciprocal geodesics
note that the requirement that its orbit be traversed an even number of times
is necessary to guarantee that the geodesic corresponds to an element of the
modular group. Namely, a lift of such a geodesic corresponds to the axis of a
hyperbolic element which passes through an order two fixed point. Since there
is only one conjugacy class of order two fixed points in $PSL(2,\mathbb{Z})$,
this hyperbolic element can be conjugated so that its axis passes through $i$,
the fixed point of $A$. As a consequence, a word representing a reciprocal
element (as an abstract group element) is conjugate to a word of the form
$a(\gamma a\gamma^{-1})$ and thus its translation length is realized on the
surface as a geodesic that retraces its steps. In terms of the generators of
$G=\mathbb{Z}_{2}\ast\mathbb{Z}_{3}$ it has the normal form,
$ab^{\epsilon_{1}}...ab^{\epsilon_{t}}ab^{-\epsilon_{t}}...ab^{-\epsilon_{1}}$,
where $\epsilon_{i}=\pm 1$. Note that the word length of a reciprocal geodesic
is always a multiple of 4. The reciprocal word has a cusp excursion of depth
bigger than $D$ if there are at least (D+1) consecutive $\epsilon_{i}$ of the
same sign. Table 1 has the association between various classes of reciprocal
geodesics, the types of compositions they correspond to, and their growth
rate.

[.5]0 Reciprocal Geodesics | Bijection To | Cardinality
---|---|---
[.5]0 Length $4t$ |  $C_{t}$: compositions of $t$ | $=2^{t-1}$
[.5]0 Length $4t$ in $S_{D}=\mathcal{E}_{0}^{D}$ |  $C_{t,D}$: Compositions of $t$ parts bounded by $D$ | $\sim\left(\frac{\alpha_{D}-1}{2+(D+1)(\alpha_{D}-2)}\right)\alpha_{D}^{t}$
[.5]0 Length $4t$ in $\mathcal{E}_{2n}^{1}$ |  $C_{t}^{2n,1}:$ $2n$ excursions of depth $>1$ | =${t}\choose{2n}$
[.5]0 Length $4t$ in $\mathcal{E}_{2}^{D}$ |  $C_{t}^{2,D}$: $2$ excursions of depth $>D$ |  $\sim\left(\frac{d_{D}^{2}}{\alpha_{D}^{D}\left(\alpha_{D}-1\right)}\right)t\alpha_{D}^{t}$

Table 1. Various classes of reciprocal geodesics and their growth

It is well-known that the primitive elements in the modular group $G$ are
comprised of

  * •

one conjugacy class of order two corresponding to the element $a$,

  * •

one conjugacy class of order three corresponding to the element $b$ and its
inverse,

  * •

one conjugacy class of parabolic elements corresponding to the element $ab$
and its inverse, and

  * •

infinitely many conjugacy classes of maximal cyclic hyperbolic subgroups.

We remind the reader that a hyperbolic element in $PSL(2,\mathbb{Z})$ that
represents a reciprocal geodesic on $S$ has axis passing through an order two
fixed point. As a result, such a hyperbolic element is conjugate to itself in
$PSL(2,\mathbb{Z})$.

Recall that $C_{t,D}$ denotes the set of compositions of $t$ with parts
bounded by $D$. Counting this set involves using the recursion relation:
$|C_{t,D}|=\sum_{i=1}^{D}|C_{t-i,D}|$ for which details can be found in [8].
We have

(2.1)
$|C_{t,D}|=rnd\bigg{(}\frac{\alpha_{D}-1}{2+(D+1)(\alpha_{D}-2)}\alpha_{D}^{t}\bigg{)}$

where $rnd(x)=\lfloor{x+\frac{1}{2}}\rfloor$ and $\alpha_{D}$ is the unique
positive root of the polynomial $z^{D}-z^{D-1}-\dots-1$ (see [8] and [2] for
details). We remark that the coefficient
$\frac{\alpha_{D}-1}{2+(D+1)(\alpha_{D}-2)}$ in equation (2.1) only depends on
$D$ and thus we denote it by $d_{D}$. We note that
$2(1-2^{-D})\leq\alpha_{D}<2$ and the $\alpha_{D}$ are increasing as $D$
increases. For the details see [8].

##  3\. Growth rates of reciprocal geodesics

Let $n$ be some non-negative integer, and $D\geq 1$ a positive integer. The
focus of this work is on reciprocal geodesics with $2n$ excursions of depth
bigger than $D$. Namely, we investigate the growth rates of such reciprocal
geodesics and how they depend on the parameters $n$ and $D$. Since reciprocal
geodesics always have an even number of excursions, we fix a positive integer
$n$ and consider reciprocal geodesics with exactly $2n$ excursions of depth
greater than $D$, and word length $4t$.

We begin by setting up a bijection between reciprocal geodesics of length $4t$
on the modular surface, and compositions of $t$. Let $\mathcal{N}_{4t}$ be the
set of normal forms for a reciprocal word of length $4t$ in
$G=\mathbb{Z}_{2}\ast\mathbb{Z}_{3}$, that is, elements of the form
$ab^{\epsilon_{1}}\dots ab^{\epsilon_{t}}ab^{-\epsilon_{t}}\dots
ab^{-\epsilon_{1}},\epsilon_{i}=\pm 1$. Let $N_{4t}$ be the set of conjugacy
classes of elements in $\mathcal{N}_{4t}$, i.e. the set of reciprocal
geodesics with length $4t$. For each conjugacy class in $N_{4t}$ there are
exactly two distinct representatives from $\mathcal{N}_{4t}$, see either
Sarnak [33] or Lemma 3.6 of [2]. Thus, the map $\xi:\mathcal{N}_{4t}\to
N_{4t}$, which sends a normal form to its conjugacy class, is 2-1. Next, let
$X_{t}=\\{(\epsilon_{1},\dots,\epsilon_{t}):\epsilon_{i}=\pm 1\\}$, i.e. the
set of binary $t$-tuples. Then the map $\Phi:\mathcal{N}_{4t}\to X_{t}$, given
by $\Phi(ab^{\epsilon_{1}}\dots ab^{\epsilon_{t}}ab^{-\epsilon_{t}}\dots
ab^{-\epsilon_{1}})=(\epsilon_{1},\dots,\epsilon_{t})$, is 1-1. Now if we
projectivize $X_{t}$, that is $PX_{t}=X_{t}/\\{\pm 1\\}$ and an element of
$PX_{t}$ is the class,
$\\{(\epsilon_{1},\dots,\epsilon_{t}),(-\epsilon_{1},\dots,-\epsilon_{t})\\}$,
we see that the map $\tau:X_{t}\to PX_{t}$ which sends a $t$-tuple to its
corresponding class, is 2-1. Lastly, the set of projectivized binary
$t$-tuples is in one to correspondence with the compositions of $t$. These
correspondences are summarized in figure 2. Thus it follows that
$|N_{4t}|=\frac{1}{2}|\mathcal{N}_{4t}|=\frac{1}{2}|X_{t}|=\frac{1}{2}(2|PX_{t}|)=|PX_{t}|=|C_{t}|$,
and we have shown that reciprocal geodesics are in 1-1 correspondence with
compositions, $C_{t}$. In fact, more is true. This 1-1 correspondence persists
if we pass to the subclasses of reciprocal geodesics which are of interest in
this article. Namely, the $D$-low lying reciprocal geodesics are in bijection
with compositions with parts bounded by $D$, and the reciprocal geodesics with
$2n$ excursions of depth bigger than $D$ (that is, $\mathcal{E}_{2n}^{D}$) are
in 1-1 correspondence with compositions with exactly $n$ parts that are bigger
than $D$ (that is, $C_{t}^{n,D}$).

${\mathcal{N}_{4t}}$${X_{t}}$${N_{4t}}$${PX_{t}}$${C_{t}}$$\scriptstyle{1-1}$$\scriptstyle{2-1}$$\scriptstyle{2-1}$$\scriptstyle{1-1}$
Figure 2. Correspondence between $N_{4t}$ and $C_{t}$

###  3.1. The case of $2n$ excursions of depth bigger than 1

Denote the subset of $N_{4t}$ consisting of reciprocal words with exactly $2n$
excursions of depth greater than one by $N_{4t}^{(2n)}$. The corresponding
subset of $PX_{t}$ consists of the classes of $t$-tuples in which there are
exactly $n$ runs of $(+1)$ or $(-1)$ of length at least two starting from the
left. Denote this subset by $PX_{t}^{(n)}$. (For example, the class
$\\{\pm(1,1,-1,1,-1,-1,-1)\\}\in PX_{7}^{(2)}$, since there are seven entries
with two runs of length at least two. The first is a run of two +1’s and the
second is a run of three -1’s.) Restricting the mappings from the diagram to
these subsets, we have that $|N_{4t}^{(2n)}|=|PX_{t}^{(n)}|$.

######  Theorem 3.1.

For any integer $n\geq 0$,

$\big{|}\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|=4t\\}\big{|}={{t}\choose{2n}}\text{
and
}\frac{\big{|}\\{\gamma\in\mathcal{E}_{2n}^{1}:|\gamma|=4t\\}\big{|}}{t^{2n}}\xrightarrow{t\to\infty}\frac{1}{2n!}$

###### Proof.

Fix integers $t>0,n\geq 0$. We have reduced the problem to computing
$|PX_{t}^{(n)}|$. Let $x\in PX_{t}^{(n)}$. Then $x$ has length $t$ with
exactly $n$ runs of +1 or -1 of length at least two. Each of the $n$ runs has
a starting position and ending position which amounts to choosing $2n$ of the
$t$ positions. Thus, any such $x$ arises by selecting $2n$ slots among the
$t$. That is, $|PX_{t}^{(n)}|={{t}\choose{2n}}$. It follows that
$|PX_{t}^{(n)}|\sim\dfrac{t^{2n}}{(2n)!},$ as $t\to\infty$ ($t\in\mathbb{N}$).

∎

###  3.2. The case of two excursions of depth bigger than $D$

We next focus on reciprocal geodesics with two excursions.

We, as usual, identify a reciprocal geodesic of length $4t$ with its normal
form in $PX_{t}^{1,D}$, the $t$-tuples with one run of length greater than $D$
and all other runs at most $D$. This projective binary sequence can in turn be
identified with compositions of $t$ with one part, which we denote by $r$,
greater than $D$, and all other parts at most $D$. We denote this set by
$C_{t}^{1,D}$. In order to count this set we decompose it in the following
way. Let $r\geq D+1$ be the length of the run in $C_{t}^{1,D}$, and let $k$ be
the position in the $t$-tuple $PX_{t}^{1,D}$ where the $r$ run starts. Note
that $1\leq k\leq t-r+1$.

\begin{overpic}[scale={.75}]{SlotsArg} \put(39.5,12.0){$\downarrow$}
\put(37.0,15.0){\text{$k^{th}$ slot}} \put(70.5,12.0){$\downarrow$}
\put(67.0,15.0){$(k+r)^{th}$ slot} \put(94.5,12.0){$\downarrow$}
\put(92.0,15.0){$t^{th}$ slot} \end{overpic} Figure 3. $2\leq k\leq t-r$

For $k$ and $r$ fixed and in the permitted ranges listed above, we have

(3.1) $|C_{t}^{1,D}|=|C_{k-1,D}||C_{t-k-r+1,D}|\text{ for }1\leq k\leq t-r+1$

This follows by counting the number of possible compositions to the left and
to the right of the $r$ run and then observing that there is only one possible
$r$ run starting in the $k^{th}$ slot. See Figure 3 for a schematic. Here we
use the convention that $|C_{0,D}|=1$.

We remark that up until this point the discussion works for $D\geq 1$. For the
rest of this section we assume $D\geq 2$. Summing over the admissible values
of $k$ and $r$, we have that

(3.2)
$\big{|}\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}\big{|}=\sum_{r=D+1}^{t}\sum_{k=1}^{t-r+1}|C_{k-1,D}||C_{t-k-r+1,D}|$

Plugging the appropriate values of expression (2.1) into (3.2) we obtain the
lower estimate for
$\big{|}\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}\big{|}$ as

(3.3)
$\sum_{r=D+1}^{t}\sum_{k=1}^{t-r+1}\left(d_{D}\alpha_{D}^{k-1}-\frac{1}{2}\right)\left(d_{D}\alpha_{D}^{t-k-r+1}-\frac{1}{2}\right)$

and the upper estimate as,

(3.4)
$\sum_{r=D+1}^{t}\sum_{k=1}^{t-r+1}\left(d_{D}\alpha_{D}^{k-1}+\frac{1}{2}\right)\left(d_{D}\alpha_{D}^{t-k-r+1}+\frac{1}{2}\right)$

The upper and lower estimates are of the form,

(3.5)
$\sum_{r}\sum_{k}d^{2}\alpha^{t-r}+a\sum_{r}\sum_{k}d\alpha^{k-1}+a\sum_{r}\sum_{k}d\alpha^{t-k-r+1}+\sum_{r}\sum_{k}\frac{1}{4}$

where $a=\pm\frac{1}{2}$ and, for ease of notation, we have set $d=d_{D}$ and
$\alpha=\alpha_{D}$.

We are interested in what happens to the terms of expression (3.5) as
$t\rightarrow\infty$.

To this end we first prove a lemma.

######  Lemma 3.2.

The sums in the following are over the admissible values of $r$ and $k$.

  1. (1)

$d^{2}\sum_{r}\sum_{k}\alpha^{t-r}\sim\frac{d^{2}}{\alpha^{D}(\alpha-1)}t\alpha^{t}$

  2. (2)

$d\sum_{r}\sum_{k}\alpha^{k-1}=O(\alpha^{t})$

  3. (3)

$d\sum_{r}\sum_{k}\alpha^{t-k-r+1}=O(\alpha^{t})$

###### Proof.

To ease notation, we use the convention $d=d_{D}$ and $\alpha=\alpha_{D}$. To
prove item (1), note that
$d^{2}\sum_{k}\alpha^{t-r}=d^{2}\alpha^{t-r}\left(t-r+1\right)$. Summing this
over the admissible $r$’s we have,

$d^{2}\sum_{r=D+1}^{t}\alpha^{t-r}\left(t-r+1\right)=$
$d^{2}\alpha^{t}\left[t\sum_{r=D+1}^{t}\alpha^{-r}-\sum_{r=D+1}^{t}(r-1)\alpha^{-r}\right]=$
$\left(d^{2}\sum_{r=D+1}^{t}\alpha^{-r}\right)t\alpha^{t}-\left(d^{2}\sum_{r=D+1}^{t}(r-1)\left(\frac{1}{\alpha}\right)^{r}\right)\alpha^{t}=$
$d^{2}\left(\frac{\alpha^{-D}-\alpha^{-t}}{\alpha-1}\right)t\alpha^{t}+O(\alpha^{t}).$

To prove item (2), we first note that a short computation yields

$\alpha^{-1}\sum_{k=1}^{t-r+1}\alpha^{k}=\frac{\alpha^{t-r+1}-1}{\alpha-1}.$

Now, summing over the admissible $r$’s we obtain,

$\sum_{r=D+1}^{t}\frac{\alpha^{t-r+1}-1}{\alpha-1}=$
$\left(\frac{1}{\alpha-1}\right)\left[\alpha^{t+1}\sum_{r=D+1}^{t}\alpha^{-r}-(t-D)\right]=O(\alpha^{t}).$

The proof of item (3) follows in a similar fashion to item (2). ∎

We are now ready to prove

######  Theorem 3.3.

For $D\geq 2$,

$\frac{\big{|}\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}\big{|}}{t\alpha_{D}^{t}}\xrightarrow{t\to\infty}\frac{d_{D}^{2}}{\alpha_{D}^{D}(\alpha_{D}-1)}.$

###### Proof.

Consider the terms of expression (3.5). Lemma 3.2 and the fact that the fourth
term grows like $t^{2}$, ensures that the first term, whose asymptotic is
$t\alpha^{t}$, dominates the expression as $t\rightarrow\infty$. Since the
first term does not depend on the constant $a$, we may conclude that the lower
and upper estimates of
${\big{|}\\{\gamma\in\mathcal{E}_{2}^{D}:|\gamma|=4t\\}\big{|}}$ have the same
asymptotic. ∎

## References

  * [1] A. Basmajian and M. Liu, Low-lying geodesics on the modular surface and necklaces. _Arxiv_
  * [2] A. Basmajian and R. Suzzi Valli, Combinatorial growth in the modular group. _Groups Geom. Dyn._ 16 (2022), no. 2, 683–703
  * [3] F. Boca, V. Paşol, A. Popa, and A. Zaharescu, Pair correlation of angles between reciprocal geodesics on the modular surface. _Algebra Number Theory_ 8 (2014), no. 4, 999–1035
  * [4] J. Bourgain and A. Kontorovich, Beyond expansion II: low-lying fundamental geodesics. _J. Eur. Math. Soc. (JEMS)_ 19 (2017), no. 5, 1331–1359
  * [5] J. Bourgain and A. Kontorovich, Beyond expansion III: Reciprocal geodesics. _Duke Math J._ 168 (2019), no. 18, 3413–3435
  * [6] P. Buser, Geometry and spectra of compact Riemann surfaces. Progress in Mathematics, 106. Birkhäuser Boston, Inc., Boston, MA 1992. xiv+454 pp.
  * [7] D. Calegari and J. Louwsma, Immersed surfaces in the modular orbifold. _Proc. Amer. Math. Soc._ 139 (2011), no. 7, 2295–2308
  * [8] G. Dresden and Z. Du, A simplified Binet formula for $k$-generalized Fibonacci numbers. _J. Integer Seq._ 17 (2014), no. 4, Article 14.4.7, 9 pp.
  * [9] V. Erlandsson, A remark on the word length in surface groups. _Trans. Amer. Math. Soc._ 372 (2019), no. 1, 441–445
  * [10] V. Erlandsson, H. Parlier, and J. Souto, Counting curves, and the stable length of currents. _J. Eur. Math. Soc. (JEMS)_ 22 (2020), no. 6, 1675–1702
  * [11] V. Erlandsson and J. Souto, Counting and equidistribution of reciprocal geodesics and dihedral groups. arXiv:2204.09956
  * [12] V. Erlandsson and J. Souto, Counting curves in hyperbolic surfaces. _Geom. Funct. Anal._ 26 (2016), no 3, 729–777
  * [13] J. Gilman, Informative words and discreteness. _Combinatorial group theory, discrete groups, and number theory_ , 147–155, Contemp. Math., 421, Amer. Math. Soc., Providence, RI, 2006
  * [14] J. Gilman, Primitive curve lengths on pairs of pants. _Infinite group theory_ , 141–155, World Sci. Publ., Hackensack, NJ, 2018
  * [15] J. Gilman and L. Keen, Cutting sequences and palindromes. _Geometry of Riemann surfaces_ , 194–216, London Math. Soc. Lecture Note Ser., 368, Cambridge Univ. Press, Cambridge, 2010
  * [16] J. Gilman and L. Keen, Enumerating palindromes and primitives in rank two free groups. _J. Algebra_ 332 (2011), 1–13
  * [17] J. Gilman and L. Keen, Lifting free subgroups of $PSL(2,\mathbb{R})$ to free groups. _Quasiconformal mappings, Riemann surfaces, and Teichmüller spaces_ , 109–122, Contemp. Math., 575, Amer. Math. Soc., Providence, RI, 2012
  * [18] J. Gilman and L. Keen, Word sequences and intersection numbers. _Complex manifolds and hyperbolic geometry_ (Guanajuato, 2001), 231–249, Contemp. Math., 311, Amer. Math. Soc., Providence, RI, 2002
  * [19] V. Guba and M. Sapir, On the conjugacy growth functions of groups. _Illinois J. Math._ 54 (2010), no. 1, 301–313
  * [20] A. Haas, Excursion and return times of a geodesic to a subset of a hyperbolic Riemann surface. Proc. Amer. Math. Soc. 141 (2013), no. 11, 3957–3967
  * [21] A. Haas, Geodesic cusp excursions and metric Diophantine approximation. Math. Res. Lett. 16 (2009), no. 1, 67–85
  * [22] A. Haas, Geodesic excursions into an embedded disc on a hyperbolic Riemann surface. Conform. Geom. Dyn. 13 (2009), 1–5
  * [23] A. Haas, The distribution of geodesic excursions into the neighborhood of a cone singularity on a hyperbolic 2-orbifold. Comment. Math. Helv. 83 (2008), no. 1, 1–20
  * [24] A. Haas, The distribution of geodesic excursions out the end of a hyperbolic orbifold and approximation with respect to a Fuchsian group. Geom. Dedicata 116 (2005), 129–155
  * [25] B. Marmolejo, Growth of Conjugacy Classes of Reciprocal Words in Triangle Groups. Ph.D. thesis, City University of New York, 2020
  * [26] M. V. Melián and D. Pestana, Geodesic excursions into cusps in finite-volume hyperbolic manifolds. Michigan Math. J. 40 (1993), no. 1, 77–93
  * [27] M. Mirzakhani, Growth of the number of simple closed geodesics on hyperbolic surfaces. _Ann. of Math. (2)_ 168 (2008), no. 1, 97–125
  * [28] R. Mor, Excursions to the cusps for geometrically finite hyperbolic orbifolds and equidistribution of closed geodesics in regular covers. Ergodic Theory Dynam. Systems 42 (2022), no. 12, 3745–3791
  * [29] P. Park, Conjugacy growth of commutators. _J. Algebra_ 526 (2019), 423–458
  * [30] M. Pollicott, Limiting distributions for geodesics excursions on the modular surface. Spectral analysis in geometry and number theory, 177–185, Contemp. Math., 484, Amer. Math. Soc., Providence, RI, 2009
  * [31] A. Randecker and G. Tiozzo, Cusp excursion in hyperbolic manifolds and singularity of harmonic measure.  J. Mod. Dyn. 17 (2021), 183–211
  * [32] I. Rivin, Growth in free groups (and other stories) – twelve years later. _Illinois J. Math._ 54 (2010), no. 1, 327–370
  * [33] P. Sarnak, Reciprocal geodesics. _Analytic number theory_ , 217–237, Clay Math. Proc., 7, Amer. Math. Soc., Providence, RI, 2007 MR 2362203
  * [34] C. Series, The modular surface and continued fractions. J. London Math. Soc. (2) 31 (1985), no. 1, 69–80
  * [35] B. Stratmann, A note on counting cuspidal excursions. Ann. Acad. Sci. Fenn. Ser. A I Math. 20 (1995), no. 2, 359–372
  * [36] D. Sullivan, Disjoint spheres, approximation by imaginary quadratic numbers, and the logarithm law for geodesics. Acta Math. 149 (1982), no. 3-4, 215–237
  * [37] C. Traina, The conjugacy problem of the modular group and the class number of real quadratic number fields. J. Numer Theory 21 (1985), no. 2, 176-184.
  * [38] M. Trin, Thurston’s compactification via geodesic currents: The case of non-compact finite area surfaces.  arXiv:2208.10763