# Correction to the photometric colors of the Gaia Data Release 2 with the
stellar color regression method

Zexi Niu  National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China  Haibo Yuan  Department of
astronomy, Beijing Normal University
19th Xinjiekouwai Street, Haidian District, Beijing, China  Jifeng Liu
National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China

###### Abstract

The second Gaia data release (DR2) delivers accurate and homogeneous
photometry data of the whole sky to an exquisite quality, reaching down to the
unprecedented milli-magnitude (mmag) level for the $G$, $G_{\rm RP}$, and
$G_{\rm BP}$ passbands. However, the presence of magnitude-dependent
systematic effects at the 10 mmag level limits its power in scientific
exploitation. In this work, using about half-million stars in common with the
LAMOST DR5, we apply the spectroscopy-based stellar color regression method to
calibrate the Gaia $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ colors. With an
unprecedented precision of about 1 mmag, systematic trends with $G$ magnitude
are revealed for both colors in great detail, reflecting changes in instrument
configurations. Color dependent trends are found for the $G_{\rm BP}-G_{\rm
RP}$ color and for stars brighter than $G$ $\sim$ 11.5 mag. The calibration is
up to 20 mmag in general and varies a few mmag/mag. A revised color-color
diagram of Gaia DR2 is given, and some applications are briefly discussed.

Astronomy data analysis, Fundamental parameters of stars, Stellar photometry

††journal: ApJ

##  1 Introduction

The second data release of the Gaia mission provides photometry data in the
$G$ band for approximately 1.7 billion sources and in the integrated $G_{\rm
BP}$ and $G_{\rm RP}$ bands for approximately 1.4 billion sources calibrated
to a consistent and homogeneous photometric system(Gaia Collaboration et al.,
2016, 2018). The $G$ magnitude is measured in the astrometric field and
extracted by a point-spread-function fitting with a typical formal uncertainty
under 1 mmag for stars of $G$ $<$ 16 mag. The $G_{\rm BP}$ and $G_{\rm RP}$
magnitudes are obtained by the $BP$ and $RP$ photometer with larger formal
uncertainties, about a few mmag for stars of $G$ $<$ 15 mag(Evans et al.,
2018; Riello et al., 2018). With such an enormous data volume and exquisite
data quality, Gaia is able to present a significant advance in all photometric
investigations.

To make full use of the mmag precision photometry yielded by the Gaia DR2,
photometric calibration to mmag precision is required. However, the systematic
effects at the 10 mmag level or the higher are shown by a number of
photometric validation tests. Arenou et al. (2018) firstly show a global
increase ($\sim$2 mmag/mag) of the $G-G_{\rm BP}$ residuals with increasing
$G$ magnitude, by subtracting fitted $G-G_{\rm BP}$ as a 3-order polynomial
function of $G_{\rm BP}-G_{\rm RP}$. They also make comparisons with the
external catalogues and find a similar global variation up to 10 mmag/mag.
Other further studies, using well-calibrated, high-quality spectral libraries,
Casagrande & VandenBerg (2018, hereafter CV18), Weiler (2018, hereafter
WEI18), and Maíz Apellániz & Weiler (2018, hereafter MAW18) conduct synthetic
photometry and compare it with the Gaia DR2 photometry. All of them detect a
similar tendency in the $G$ magnitude and propose a linear correction for 3.4,
3.5, and 3.2 mmag/mag separately. Moreover, WEI18 and MAW18 find a significant
color-dependent offset of 20 mmag between faint and bright stars in $G_{\rm
BP}$. Nevertheless, since only hundreds of stars are available in their
method, their linear corrections have the exact tendency yet lack of fine
details.

Taking advantage of the fact that millions of stellar spectra and their
associated precise stellar atmospheric parameters are now available, Yuan et
al. (2015c) has proposed a spectroscopy-based stellar color regression (SCR)
method to achieve mmag precision color calibration of wide-field imaging
surveys. The method was applied to the Sloan Digital Sky Survey (SDSS; York et
al. 2000) Stripe 82 data, achieving a precision $\sim$ 2 – 5 mmag for
different SDSS colors. It is straightforward and able to capture subtle
structures thanks to the intensive database. Analogous to the SDSS database,
the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Zhao et
al. 2012; Deng et al. 2012; Liu et al. 2014) spectroscopic sky survey is also
appropriate to do so. It begins in 2012, acquiring thousands of stellar
spectral per night with the spectral resolution R $\sim$ 1800 and accumulating
more than 8 million stellar spectra in DR5. Effective temperature $T_{\rm
eff}$, surface gravity log $g$, and metallicity [Fe/H], are delivered from the
spectra through the LAMOST Stellar Parameter pipeline (LASP; Luo et al., 2015;
Wu et al., 2011) with the precision of about 110 K, 0.2 dex, and 0.1 dex,
respectively (Luo et al., 2015). In this work, combing the LAMOST DR5
spectroscopy and the Gaia DR2 photometry, we apply the SCR method and reach
new corrections to the Gaia photometry at an unprecedented precision of about
1 mmag.

The paper is organized as follows. We describe our data selection in Section
2. The method and result are demonstrated in Section 3. The result validation
and the comparison with previous researches are discussed in Section 4. We
summarize our paper in Section 5.

##  2 Data

Figure 1: Top: Gaia DR2 observed $G_{\rm BP}-G_{\rm RP}$ color against excess
factor for three different $G$ magnitude bins at 1 mag interval. Five hundred
randomly selected stars are plotted for each bin. Bottom: same plot but for 3
metallicity bins of the same $G$ magnitude range 11-12 mag. The black dots are
dropped suspicious objects in the 2$\sigma$-clipping process.

Stars having high-precision photometric colors and well-determined stellar
atmospheric parameters are required in our study. Here we consider the
following constraints to guarantee the data quality:

  * 1)

phot$\\_$bp$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 2)

phot$\\_$g$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 3)

phot$\\_$rp$\\_$mean$\\_$flux$\\_$over$\\_$error $>$ 100

  * 4)

phot$\\_$proc$\\_$mode = 0

  * 5)

duplicated$\\_$source = False

  * 6)

phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor (see text for limit)

  * 7)

$E(B–V)_{\rm SFD}<$ 0.05 mag

  * 8)

$|\texttt{glat}|$ $>$ 20 $\deg$

  * 9)

$|Z|\ >$ 0.2 kpc, where Z is vertical distance to the Galactic disk

  * 10)

the signal-to-noise ratios for the g band ($S/N_{\rm g}$) of the LAMOST
spectra are larger than 30

  * 11)

$T_{\rm eff}>$ 4500 K for dwarf stars

  * 12)

remove variable stars

Among them, criteria (1)-(6) are applied to the Gaia DR2 data.
phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor is defined as the excess of flux in
the $G_{\rm BP}$ and $G_{\rm RP}$ with the respect to the $G$ band, which is
believed to be caused by background and contamination issues affecting the
$G_{\rm BP}$ and $G_{\rm RP}$ data (Evans et al., 2018). We find that it has a
significant systematic trend with apparent $G$ magnitude and [Fe/H]. As shown
in the Figure 1, the phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor systematically
increases with the decreasing $G$ magnitude (about 0.005 per magnitude) and
decreases with increasing metallicity (about 0.001 per 0.5 dex) within a given
$G$ magnitude range. We thus bin our sample in the $G$ magnitude and [Fe/H]
space with 1 mag and 0.5 dex intervals respectively and apply a
2$\sigma$-clipping on the phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor to remove
outliers in each bin (see Figure 1). About 4$\%$ of objects are excluded in
this step. Note that due to the relatively large widths of the magnitude and
[Fe/H] bins, the estimated $\sigma$ values are larger than the scatters caused
by the errors of the phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor. Therefore, a
2$\sigma$-clipping is used, a little bit excessive but secure.

We use the Schlegel, Finkbeiner, & Davis (1998, hereafter SFD) dust reddening
map and criteria (7)-(9) to select low extinction stars for reliable reddening
correction.

The empirically determined, temperature and reddening dependent reddening
coefficients R($G-G_{\rm RP}$) and R($G_{\rm BP}-G_{\rm RP}$) are given by the
following functions (Sun et al., to be submitted):

$\displaystyle R(G_{\rm BP}-G)=0.901-0.225\times E(B-V)_{\rm SFD}+0.178$ (1)
$\displaystyle\times{E(B-V)}^{2}_{\rm SFD}-6.422\times 10^{-5}\times T_{\rm
eff}-3.982\times 10^{-7}$ $\displaystyle\times T_{\rm eff}\times E(B-V)_{\rm
SFD}+3.541\times 10^{-9}\times{T_{\rm eff}}^{2}$

$\displaystyle R(G_{\rm BP}-G_{\rm RP})=1.060-0.253\times E(B-V)_{\rm
SFD}+0.264$ (2) $\displaystyle\times{E(B-V)}^{2}_{\rm SFD}+5.535\times
10^{-6}\times T_{\rm eff}-4.505\times 10^{-5}$ $\displaystyle\times T_{\rm
eff}\times E(B-V)_{\rm SFD}+6.64\times 10^{-9}\times{T_{\rm eff}}^{2}$

Figure 2: $T_{\rm eff}$ against log $g$ of sources satisfying the filters
described in Sect.2. Below the cyan line are dwarf stars from the LAMOST DR5,
while the giant stars are from the value-added RGB catalog of LAMOST DR4.

Note that all colors referred to hereafter are dereddened intrinsic ones.

Criteria (10) and (11) aim to select high-quality stellar parameters from the
LAMOST. For stars with multiple observations, their mean values using inverse
variance weighting are adopted.

For criterion (12), given that the Gaia satellite scans the sky repeatedly and
publishes a mean magnitude for each object, variable stars could be less
reliable and should be removed. A number of 112 objects identified as variable
stars in the WISE (Chen et al., 2018b), ATLAs (Heinze et al., 2018), ASAS-SN
(Jayasinghe et al., 2020), Catalina (Drake et al., 2017), and ZTF (Chen et
al., 2020) surveys are removed in this step.

Finally, we divide the above sample into dwarfs and giants, considering that
they probably suffer different systematic errors in the effective temperature
and metallicity. Dwarf stars are selected from $T_{\rm eff}$ against log $g$
diagram as those below the cyan line in Figure 2. For giants, only Red Giant
Brunch (RGB) stars are selected (with red clump stars removed) by cross-
matching the above sample with a value-added RGB catalog of LAMOST DR4 (Wu et
al., 2019). The final sample contains 433,845 dwarfs and 55,395 RGB stars, as
displayed in Figure 2.

##  3 Method and Result

Figure 3: Top: Differences of LAMOST DR5 stellar parameters deduced from two
observations of very different $S/N_{\rm g}$ as a function of the lower
$S/N_{\rm g}$. The spectrum of higher quality has $S/N_{\rm g}>100$. The red
lines indicate the medians and standard deviations. Note the very small
standard deviation values at $S/N_{\rm g}>30$. Bottom: Differences between
LAMOST DR5 and APOGEE DR16 stellar parameters for solar-type stars as a
function of the LAMOST $S/N_{\rm g}$. The APOGEE stars have spectral S/Ns
higher than 100.

The key idea of the SCR method is that stars of the same stellar atmospheric
parameters should have the same intrinsic colors, therefore, plenty of stars
with precisely determined spectroscopic parameters can serve as excellent
color standards to carry out very precise color calibrations. Yuan et al.
(2015c) demonstrates the method with the SDSS Stripe 82 data. In their work, a
well-calibrated control sample within a specific area is selected as a
reference to define the intrinsic colors as a function of the stellar
parameters, which in turn are used to map out the spatial variations of the
color zero-point offset as functions of R.A. and Dec., achieving a remarkable
improvement in color calibration by a factor of two to three. Very recently,
Huang et al. (2020) apply the SCR method to the second data release of the
SkyMapper Southern Survey (Onken et al., 2019). A uniform calibration with
precision better than 1$\%$ is achieved.

In this work, taking advantages of about half million common stars selected in
the previous section, we apply this method to the Gaia DR2 and obtain the
empirical relations between the Gaia colors ($G-G_{\rm RP}$ and $G_{\rm
BP}-G_{\rm RP}$) and the LAMOST DR5 stellar astrophysical parameters ($T_{\rm
eff}$, [Fe/H], log $g$). Since a systematic effect with the $G$ magnitude is
known, a control sample within a narrow magnitude range is needed as the
reference to explore the magnitude-dependent variations. Gaia modifies its
instrument configurations (Gates and Windows) according to the $G$ magnitude,
even though an iterative calibration is performed to establish the internal
photometric system, the convergence problems still exist among the different
configurations, and stars around these joints may have relatively larger
uncertainties (Evans et al., 2018). Considering the distribution of the Gates
and Window Classes (Evans et al., 2018) as well as the number of stars, stars
of $13.3<G<13.7$ are selected as the control sample, consisting of 50,343
dwarfs and 6,144 RGB stars.

Figure 4: Fitting residuals of the $G-G_{\rm RP}$ (top) and $G_{\rm BP}-G_{\rm
RP}$ (bottom) of the control sample for the dwarf stars as functions of
$T_{\rm eff}$, [Fe/H], and log $g$. The red lines show the medians and
standard deviations of the residuals. Figure 5: Same as Figure 4, but for the
RGB stars.

Magnitude-dependent systematic effects in the LAMOST stellar parameters could
cause false trends with $G$ magnitudes in our results. We select multiply-
observed stars in the LAMOST DR5 to investigate such effects and find that
there are no systematic effects with the $S/N_{\rm g}$ till at very low
$S/N_{\rm g}$ of about 15. The results are plotted in the top panels of Figure
3 showing a comparison between observations of very different $S/N_{\rm g}$.
To further verify this result, we cross-match the LAMOST DR5 with the APOGEE
DR16 (Ahumada et al., 2020) and select a sample of solar-type stars to plot
their stellar parameter differences as a function of LAMOST $S/N_{\rm g}$ in
the bottom panels of Figure 3. A flat offset is seen for each stellar
parameter at $S/N_{\rm g}>15$, consistent with the test result from duplicated
stars. Note the offsets are not exactly zero due to systematic errors in the
LAMOST and APOGEE stellar parameters.

To precisely construct the relationships between the Gaia colors and the
LAMOST stellar parameters, we divide our control sample into different bins as
following: $T_{\rm eff}$ in 300 K, [Fe/H] in 0.2 dex, and log $g$ in 0.3 dex
for dwarfs; $T_{\rm eff}$ in 300 K, [Fe/H] in 0.5 dex, and log $g$ in 0.4 dex
for RGB stars. A linear fitting formula is adopted for each grid:
$C=a_{0}+a_{1}\times T_{\rm eff}+a_{2}\times[Fe/H]+a_{3}\times$ log $g$, where
C can be ($G_{\rm BP}-G_{\rm RP}$)0 or ($G-G_{\rm RP}$)0. 2$\sigma$-clipping
is performed during the fitting process. Bins containing less than 50 stars
are excluded to avoid poor fitting. The fitting residuals are plotted in
Figure 4 and Figure 5. The median residuals of both dwarfs and RGB giants are
very close to zero. The typical standard deviations are about 7.5 mmag and 15
mmag for $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$, respectively.

Figure 6: Distributions of the color residuals (fitted colors $-$ observed
colors) in the $T_{\rm eff}$ and ($G_{\rm BP}-G_{\rm RP}$)0 plane for the
dwarf stars. The black solid lines mark the boundaries of different
subsamples. To avoid the selection effect, only colorful points are used. Only
1 in 100 stars are plotted.  Figure 7: Same as Figure 6, but for the RGB
stars. Only 1 in 10 stars are plotted.

Applying the empirical relations derived above to the whole sample, color
residuals/corrections (fitted colors $-$ observed colors) are obtained. Their
distributions in the $T_{\rm eff}$ and ($G_{\rm BP}-G_{\rm RP}$)0 plane are
shown in Figure 6 and Figure 7. At a given effective temperature, stars of
redder intrinsic colors have typically smaller (negative) residuals. To
examine the color dependence, both dwarfs and RGB giants are divided into
subsamples by the black solid lines in Figure 6 and Figure 7. Four subsamples
are adopted for dwarfs (Main sequence, hereafter MS), with typical ($G_{\rm
BP}-G_{\rm RP}$)0 colors of 1.06, 0.90, 0.77, and 0.66 mag, respectively. Two
subsamples are adopted for RGB stars, with typical ($G_{\rm BP}-G_{\rm RP}$)0
colors of 1.13 and 1.01 mag, respectively. The sample numbers from red to blue
are 32655, 102441, 163325, 68493 for dwarfs and 24697, 14601 for giants. The
gray points are discarded because stars there are partly dropped due to the
requirement N $>$ 50, and the dramatic change of numbers at the margin can
cause bias of color residuals.

For each subsample, the distributions of its color residuals against $G$
magnitudes are plotted in Figure 8 and Figure 9. The points are grouped into
bins along the $G$ magnitude for every 0.05 mag. The median and standard
deviation values are over-plotted in orange lines. Typical standard deviations
are $\sim$ 7 – 9 mmag for $G-G_{\rm RP}$ and $\sim$ 12 – 17 mmag for $G_{\rm
BP}-G_{\rm RP}$. Magnitude dependent offsets are clearly seen for all
subsamples. The offsets (median of the color residuals) are smoothed by the
locally estimated scatterplot smoothing (LOWESS) and taken as the calibration
curves. The algorithm works by taking the $frac\times N$ closest points to
each data point and estimating smoothed $y$ values using a weighted linear
regression based on their $x$ distances, where $N$ is the total number of
points and $frac$ varies for subsamples with typically 0.05. Note that the
curves of different subsamples may have different magnitude ranges. In places
with few sources at the very bright or faint ends, linear extrapolation is
performed. Specifically, corrections can be expressed as:

$C^{{}^{\prime}}=C+\Delta C$ (3)

where $C^{{}^{\prime}}$ is the corrected color, $C$ is the Gaia DR2 color and
$\Delta C$ is the color correction term.

Figure 8: Fitting residuals of $G-G_{\rm RP}$ against $G$ magnitude for each
subsample from Figure 6 and Figure 7 , labeled by its evolution stage and
typical ($G_{\rm BP}-G_{\rm RP}$)0 color. Each subsample is binned along the
$G$ magnitude with 0.05 mag interval. The median and standard deviation values
are estimated and shown by the orange lines. A 2$\sigma$-clipping is performed
in this process. The changes of instrument configurations are marked by the
vertical blue lines.

In Figure 8, there are three marked variations at around $G$ $\sim$ 11.2,
13.0, and 15.7 mag that are probably related to the changes of instrument
configurations. Gates are implemented on the Gaia CCDs to ease the saturation
issues in the case of bright stars, e.g. $G$ $<$ 12 mag. Sizes of Windows in
along-scan and across-scan direction are modified at $G$ = 13.0 and 16 mag for
the $G$ band, and at $G$ = 11.5 for the $G_{\rm BP}$ and $G_{\rm RP}$ bands
(Evans et al., 2018). However, the choice of the instrument configurations
depends on the observed signals in the $G$ band, the corresponding calibrated
$G$ magnitudes are not exactly identical for every single observation. For
example, the transition of the $G$ band from Window Class 0 to 1 may move
around $G$ = 13 mag, and the transition of the $G_{\rm BP}$ and $G_{\rm RP}$
may move around $G$ = 11.5 mag, due to issues including flat-fielding and
stray light contamination. Since the homogeneous calibration adopted by Gaia
is a function of the $G$ magnitude and may differ significantly before and
after the transition magnitude, some corrections could be incorrect.
Consequently, the random uncertainties around the transitions are larger than
those in adjacent regions, as demonstrated by Figure 9 in Evans et al. (2018).
The features at $G$ $\sim$ 11.2, 13.0, and 15.7 mag in Figure 8 are consistent
with Figure 9 in Evans et al. (2018), suggesting that these fine structures
are real and accurate.

Figure 9: Same as Figure 9, but for the $G_{\rm BP}-G_{\rm RP}$ color.

The final results and a comparison of calibrations of different subsamples are
plotted in Figure 10. Positions of changes of Gates and Windows are marked as
well. The analysis and possible applications of our results are discussed in
Section 4.

Figure 10: Color-dependent and recommended calibration curves for $G-G_{\rm
RP}$ and $G_{\rm BP}-G_{\rm RP}$ .  Figure 11: Comparisons of spatial
distributions in Mollweide projection of the dwarfs and RGB giants of 12 $<G<$
13 mag (top) and 15 $<G<$ 16 mag (bottom).

##  4 Discussion

###  4.1 Color dependence

As shown in Figure 10, despite places with few sources and relatively larger
errors at the very bright ends, the calibration curves of $G-G_{\rm RP}$
yielded by the six different subsamples agree well with each other when $G$
$<$ 14 mag. However, when $G$ $>$ 14 mag, discrepancies exist between the
curves of the two RGB subsamples (RGB 1.13 and RGB 1.01) and the red dwarf
subsample (MS 1.06), while those of the other three blue dwarf subsamples (MS
0.90, MS 0.77, and MS 0.66) still agree well. As for the calibration curves
for $G_{\rm BP}-G_{\rm RP}$, the same trend exists as well when $G$ $>$ 14
mag. When $G$ $<$ 11.5 mag, the curves are divided into two groups (RGB 1.13
and RGB 1.01 versus MS 0.77 and MS 0.66) according to their colors.

The discrepancies of the three red lines (MS 1.06, RGB 1.13, and RGB 1.01) at
$G$ $>$ 14 mag in Figure 10 are not real, but caused by the combination of two
reasons: the selection function of the LAMOST data and the spatially dependent
systematics of the SFD reddening map. The LAMOST surveys have four different
types of plates: very bright (VB), bright (B), median bright (M), and faint
(F) plates, in which VB plates target stars of $r<$ 14 mag while the others
target fainter stars (Chen et al., 2018a). Compared to blue stars, bright red
stars are more likely to be giants while faint red stars to be dwarfs. The
comparisons of the spatial distributions of dwarfs and RGB stars are shown in
Figure 11. The whole sky is divided into different bins of 10 degree by 10
degree in size. The percentages $f_{i}=\frac{N_{i}}{N_{\rm total}}$ of RGB
stars and dwarfs in each bin are calculated ,and the colorbars represent the
ratios of dwarf percentage to RGB percentage $\frac{f_{i\rm RGB}}{f_{i\rm
MS}}$. The closer the ratios get to 1, the smaller the differences between the
spatial distributions of RGB stars and dwarfs are. Clearly, the differences
are larger when $G>$ 14 mag. Note that the Gaia $G$ magnitudes are very close
to the SDSS $r$ magnitudes for typical stars.

On the other hand, the SFD reddening map is found to show systematic errors
that depend on spatial position and dust temperature (e.g., Peek & Graves,
2010; Sun et al., to be submitted). The top and middle panels of Figure 12
show sky distributions of $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ residuals
(fitted colors $-$ observed colors) after correcting for the trends with $G$
magnitude, which are almost identical and therefore caused by the same
problem. The SFD map overestimates reddening in the pink regions (positive
residuals) while underestimates reddening in the blue regions (negative
residuals), which is consistent with the result of Sun et al. (to be
submitted) well. A strong correlation is seen in the the bottom panel of
Figure 12 that compares the corrected $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm
RP}$ residuals of the whole sample. The red crosses are the median values of
two corrected residuals. They agree very well with the cyan line, whose slope
is determined by the median value ($\sim$ 1.99) of the $\frac{R(G_{\rm
BP}-G_{\rm RP})}{R(G-G_{\rm RP})}$ ratios of the sample, as expected.

For the $G_{\rm BP}$ and $G_{\rm RP}$ bands, even though the Gaia photometric
pipeline treats both 2D and 1D Windows as aperture photometry in the same way
(Evans et al., 2018), different instrument configurations may still cause
problems. WEI18 points out that convergence during the calibration process may
be dominated by stars having the most number, which can cause discrepancies
among colors. They find that $G_{\rm BP}$ magnitude has a systematic tendency
among $G_{\rm BP}-G_{\rm RP}$ and their offset happens around 10.99 mag. MAW18
further confirms that the color-dependent break happens in $G_{\rm BP}$ band
at $G$ $\sim$ 10.87 mag, by analyzing the uncertainties of the $G_{\rm BP}$
fluxes. Likewise, we believe that the color-dependent feature at $G$ $<$ 11.5
mag in our result of $G_{\rm BP}-G_{\rm RP}$ calibration is due to the same
reason.

We combine the three unbiased dwarf subsamples (MS 0.90, MS 0.77, and MS 0.66)
together to determine the recommended calibration curves with the same
procedure in Section 3. The recommended calibration curves are overplotted in
black lines in Figure 10. The differences between the recommended curve and
the three colored ones (MS 0.90, MS 0.77, and MS 0.66) are computed, with
standard deviations 0.3 and 0.6 mmag for $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm
RP}$, respectively. This suggests that the calibration curve of each subsample
has a typical random error smaller than 1.0 mmag. Note that the recommended
curves work well in most cases except for $G_{\rm BP}-G_{\rm RP}$ when $G$ $<$
11.5 mag, which is plotted in the dash-dot line. The recommended color
calibration curves as well as the calibration curves of different subsamples
are listed in Table 2.

Table 1: Color calibration curves. The first column is $G$ magnitude. The second and third columns are respectively recommended $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ calibration curves, followed by those yielded by different subsamples. All calibration curves are in units of mmag.  G | recom | MS 1.06 | MS 0.89 | MS 0.77 | MS 0.67 | RGB 1.13 | RGB 1.02
---|---|---|---|---|---|---|---
… | … | … | … | … | … | … | … | … | … | … | … | … | … | …
13.01 | 0.57 | 1.53 | 0.40 | 0.43 | 0.34 | 1.37 | 0.53 | 1.46 | 1.07 | 1.84 | 1.13 | 1.27 | 1.52 | 2.05
13.02 | 0.46 | 1.50 | 0.40 | 0.46 | 0.28 | 1.34 | 0.46 | 1.42 | 0.94 | 1.78 | 1.06 | 1.24 | 1.44 | 2.00
13.03 | 0.37 | 1.46 | 0.40 | 0.49 | 0.22 | 1.31 | 0.39 | 1.38 | 0.81 | 1.72 | 0.99 | 1.21 | 1.36 | 1.95
13.04 | 0.28 | 1.43 | 0.40 | 0.52 | 0.16 | 1.28 | 0.32 | 1.34 | 0.68 | 1.66 | 0.92 | 1.18 | 1.28 | 1.90
13.05 | 0.21 | 1.39 | 0.40 | 0.55 | 0.10 | 1.25 | 0.25 | 1.30 | 0.55 | 1.60 | 0.85 | 1.15 | 1.20 | 1.85
13.06 | 0.14 | 1.36 | 0.40 | 0.58 | 0.04 | 1.22 | 0.18 | 1.26 | 0.42 | 1.54 | 0.78 | 1.12 | 1.12 | 1.80
13.07 | 0.09 | 1.33 | 0.40 | 0.61 | -0.02 | 1.19 | 0.11 | 1.22 | 0.29 | 1.48 | 0.71 | 1.09 | 1.04 | 1.75
13.08 | 0.05 | 1.29 | 0.40 | 0.64 | -0.08 | 1.16 | 0.04 | 1.18 | 0.16 | 1.42 | 0.64 | 1.06 | 0.96 | 1.70
13.09 | 0.01 | 1.26 | 0.40 | 0.67 | -0.14 | 1.13 | -0.03 | 1.14 | 0.03 | 1.36 | 0.57 | 1.03 | 0.88 | 1.65
13.10 | -0.01 | 1.23 | 0.40 | 0.70 | -0.20 | 1.10 | -0.10 | 1.10 | -0.10 | 1.30 | 0.50 | 1.00 | 0.80 | 1.60
13.11 | -0.02 | 1.19 | 0.42 | 0.74 | -0.19 | 1.07 | -0.09 | 1.07 | -0.08 | 1.26 | 0.47 | 0.98 | 0.76 | 1.57
13.12 | -0.02 | 1.16 | 0.44 | 0.78 | -0.18 | 1.04 | -0.08 | 1.04 | -0.06 | 1.22 | 0.44 | 0.96 | 0.72 | 1.54
13.13 | -0.02 | 1.13 | 0.46 | 0.82 | -0.17 | 1.01 | -0.07 | 1.01 | -0.04 | 1.18 | 0.41 | 0.94 | 0.68 | 1.51
13.14 | -0.01 | 1.09 | 0.48 | 0.86 | -0.16 | 0.98 | -0.06 | 0.98 | -0.02 | 1.14 | 0.38 | 0.92 | 0.64 | 1.48
13.15 | 0.00 | 1.06 | 0.50 | 0.90 | -0.15 | 0.95 | -0.05 | 0.95 | -0.00 | 1.10 | 0.35 | 0.90 | 0.60 | 1.45
13.16 | 0.02 | 1.03 | 0.52 | 0.94 | -0.14 | 0.92 | -0.04 | 0.92 | 0.02 | 1.06 | 0.32 | 0.88 | 0.56 | 1.42
13.17 | 0.03 | 1.00 | 0.54 | 0.98 | -0.13 | 0.89 | -0.03 | 0.89 | 0.04 | 1.02 | 0.29 | 0.86 | 0.52 | 1.39
13.18 | 0.04 | 0.96 | 0.56 | 1.02 | -0.12 | 0.86 | -0.02 | 0.86 | 0.06 | 0.98 | 0.26 | 0.84 | 0.48 | 1.36
13.19 | 0.05 | 0.93 | 0.58 | 1.06 | -0.11 | 0.83 | -0.01 | 0.83 | 0.08 | 0.94 | 0.23 | 0.82 | 0.44 | 1.33
13.20 | 0.06 | 0.90 | 0.60 | 1.10 | -0.10 | 0.80 | -0.00 | 0.80 | 0.10 | 0.90 | 0.20 | 0.80 | 0.40 | 1.30
… | … | … | … | … | … | … | … | … | … | … | … | … | … | …

###  4.2 Comparisons with previous works

Evans et al. (2018) publish a set of revised sensitivity curves called REV,
which is preferred but have no significant differences with the DR2 ones. In
CV18, with the well-calibrated CALSPEC spectral library, they compute the
synthetic $G$, $G_{\rm BP}$, and $G_{\rm RP}$ magnitudes using the REV curves
and compare with those in the Gaia DR2 catalog, yielding a linear correction
of 3.5 mmag/mag in the $G$ magnitude. WEI18 reconstructs the Gaia passbands in
basic functional analytic framework with four spectral libraries including the
CALSPEC and finds a close magnitude drift with CV18. Later, MAW18 updates the
spectral libraries and applies the same method in WEI18, suggesting a 3.2
mmag/mag correction in the $G$ magnitude.

Figure 13 compares our results with those of CV18 and MAW18. For the $G-G_{\rm
RP}$ color, systematic offsets can be seen between any two of the three
references because of the different processing progress: (1) we take $G$ =
13.5 $\pm$ 0.2 mag as the control sample; (2) CV18 applies the REV
transmission curves and the official zero-points; and (3) MAW18 uses their own
transmission curves called MAW and a zero-point of 0. Both CV18 and MAW18
provide a linear fitting of the systematic trend in the $G$ band. Our result
shows a very similar trend, consistent with the results of CV18 and MAW18.
However, thanks to a much larger sample of stars used (half million compared
to about one hundred), more fine details are revealed in this work,
particularly those due to the changes of the instrument configurations as
discussed in Section 3. For the $G_{\rm BP}-G_{\rm RP}$ color, our result is
also consistent with those of CV18 and MAW18 within their errors. Our result
shows a small but clear trend with G magnitude at G $>$ 11.5 mag. At the
bright end, our result is slightly higher, probably due to the fact that most
stars there in CV18 and MAW18 are bluer.

Figure 12: Bias of the SFD dust reddening map. Top: Scatter plot of the sky
distribution in Mollweide projection of the corrected $G-G_{\rm RP}$ residual
(fitted colors $-$ observed colors). Middle: Same but for $G_{\rm BP}-G_{\rm
RP}$. Only 10000 randomly selected points are plotted. Bottom: Comparison of
the corrected color residuals. The red crosses are the median values of two
corrected residuals. Instead of fitting the red crosses, the slope of the cyan
line is determined by the median value of the $\frac{R(G_{\rm BP}-G_{\rm
RP})}{R(G-G_{\rm RP})}$ ratios of the sample stars.  Figure 13: Comparisons of
our results (red lines) with those of MAW18 (empty circles) and CV18 (black
crosses). Note that the typical errors are about 10 mmag for the MAW18 and
CV18 data.

###  4.3 Revised color-color diagram

Here we carry out a test of the validity of our calibration curves by
comparing the observed and the revised color–color diagrams. The recommended
blank solid curves are adopted to correct $G-G_{\rm RP}$ and $G_{\rm
BP}-G_{\rm RP}$ for stars of $G$ $>$ 11.5 mag. For stars of $G$ $<$ 11.5 mag,
$G_{\rm BP}-G_{\rm RP}$ colors are corrected by colored curves. For stars of
$G_{\rm BP}-G_{\rm RP}$ $>$ 1.13 mag, the calibration curve of $G_{\rm
BP}-G_{\rm RP}$ = 1.13 mag is used. We correct the whole sample using equation
(3) and compare color-color diagrams before and after correction in the left
and middle panels in Figure 14. To exhibit the metallicity-dependent stellar
locus clearly, only 300 stars in each metallicity bin are plotted and colored
by their metallicities. After correction, the distribution of stars with
similar [Fe/H] gets tighter, as expected. Following the work of Yuan et al.
(2015a), a fourth-order polynomial is adopted to fit the metallicity-dependent
stellar loci of the Gaia colors. The histograms of the fitting residuals are
plotted in the right panels of Figure 14. Both dwarfs and giants show a much
smaller $\sigma$ after correction, from 2.51 mmag to 1.76 mmag for dwarfs and
from 2.59 mmag to 1.55 mmag for giants. The distribution for giants also
becomes more symmetric and Gaussian.

Like Figure 31 in Arenou et al. (2018), we also plot a 2D histogram of the
$G-G_{\rm BP}$ residuals in Figure 15 after subtracting the metallicity-
dependent stellar color locus. Subplot labelled DR2 is identical to Arenou et
al. (2018), showing a clear trend with G magnitude. After the linear
correction in the G magnitude proposed by CV18, the overall trend disappears,
but still shows small scale features at a few mmag level. Likewise, the liner
corrections suggested by WEI18 and MAW18 suffer the same problem. After
corrections of this work, the trend is flat everywhere.

The revised color-color diagram is essential in a number of studies, for
instance, detecting peculiar objects or determining reliable metallicities for
an enormous and magnitude-limited sample of stars from the Gaia photometry (Xu
et al., to be submitted). We notice that the $G_{\rm BP}-G$ residual
distribution of dwarfs after correction is asymmetric at its wings. There is
an excess of stars whose colors are redder compared to those predicted by the
metallicity-dependent color locus. This asymmetric features can be well
explained by binary stars (Yuan et al., 2015b) and used to estimate binary
fractions of a huge number of field stars (Niu et al., to be submitted).

Figure 14: Left and middle columns compare $G_{\rm BP}-G_{\rm RP}$ vs. $G_{\rm
BP}-G$ diagrams before and after corrections, respectively. Right column
compares the $G_{\rm BP}-G$ residual distributions before and after
corrections. Gaussian fitting is applied, and the fitted $\sigma$ values are
texted. top: dwarfs; bottom: RGB giants. Figure 15: 2D histogram of the
$G-G_{\rm BP}$ residual after subtracting the metallicity-dependent color
locus. top: published DR2 data. middle: applying CV18 corrections. bottom:
applying corrections from this work.

##  5 Summary

In this work, using about half-million stars with high-quality Gaia DR2
photometry and LAMOST DR5 stellar parameters , we apply the SCR method to
calibrate the Gaia $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ colors. An
unprecedented precision of about 1 mmag is achieved, suggesting the great
power of the SCR method in calibrating photometric surveys.

Magnitude dependent trends are revealed for both the $G-G_{\rm RP}$ and
$G_{\rm BP}-G_{\rm RP}$ colors in great details, reflecting changes in
instrument configurations. Color dependent trends are found for the $G_{\rm
BP}-G_{\rm RP}$ color and for stars brighter than $G$ $\sim$ 11.5 mag. The
calibration is up to 20 mmag in general and varies a few mmag/mag. Our results
are consistent with previous results but with a much improved precision.
Implementations of the color and magnitude corrections are provided. A revised
color-color diagram of Gaia DR2 is given to demonstrate some potential
applications of the corrected mmag photometry. The $G$, $G_{\rm BP}$, and
$G_{\rm RP}$ passbands are different in Gaia EDR3 in the same way that the $G$
passbands were different between DR1 and DR2. We will extend the work to Gaia
EDR3 data in a separate paper.

We acknowledge the anonymous referee for his/her valuable comments that
improve the quality of this paper significantly. We acknowledge Profs. L.
Casagrande and J. Maíz Apellániz for providing data of their work. We greatly
thank Huiqin Yang for a careful reading of the manuscript. This work is
supported by National Science Foundation of China (NSFC) under grant numbers
11603002, 11988101, and 113300034, National Key Research and Development
Program of China (NKRDPC) under grant numbers 2016YFA0400804, 2019YFA0405503,
and 2019YFA0405504. This work has made use of data products from the
Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope, LAMOST). LAMOST is a National Major Scientific Project built by the
Chinese Academy of Sciences. Funding for the project has been provided by the
National Development and Reform Commission. LAMOST is operated and managed by
the National Astronomical Observatories, Chinese Academy of Sciences. This
work has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Ahumada et al. (2020) Ahumada, R., Allende Prieto, C., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Arenou et al. (2018) Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17, doi: 10.1051/0004-6361/201833234
  * Casagrande & VandenBerg (2018) Casagrande, L., & VandenBerg, D. A. 2018, MNRAS, 479, L102, doi: 10.1093/mnrasl/sly104
  * Chen et al. (2018a) Chen, B. Q., Liu, X. W., Yuan, H. B., et al. 2018a, MNRAS, 476, 3278, doi: 10.1093/mnras/sty454
  * Chen et al. (2018b) Chen, X., Wang, S., Deng, L., de Grijs, R., & Yang, M. 2018b, ApJS, 237, 28, doi: 10.3847/1538-4365/aad32b
  * Chen et al. (2020) Chen, X., Wang, S., Deng, L., et al. 2020, ApJS, 249, 18, doi: 10.3847/1538-4365/ab9cae
  * Deng et al. (2012) Deng, L.-C., Newberg, H. J., Liu, C., et al. 2012, Research in Astronomy and Astrophysics, 12, 735, doi: 10.1088/1674-4527/12/7/003
  * Drake et al. (2017) Drake, A. J., Djorgovski, S. G., Catelan, M., et al. 2017, MNRAS, 469, 3688, doi: 10.1093/mnras/stx1085
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4, doi: 10.1051/0004-6361/201832756
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Heinze et al. (2018) Heinze, A. N., Tonry, J. L., Denneau, L., et al. 2018, AJ, 156, 241, doi: 10.3847/1538-3881/aae47f
  * Huang et al. (2020) Huang, Y., Yuan, H., Li, C., et al. 2020, arXiv e-prints, arXiv:2011.07172.  https://arxiv.org/abs/2011.07172
  * Jayasinghe et al. (2020) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2020, MNRAS, 491, 13, doi: 10.1093/mnras/stz2711
  * Liu et al. (2014) Liu, X. W., Yuan, H. B., Huo, Z. Y., et al. 2014, in Setting the scene for Gaia and LAMOST, ed. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, Vol. 298, 310–321, doi: 10.1017/S1743921313006510
  * Luo et al. (2015) Luo, A. L., Zhao, Y.-H., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1095, doi: 10.1088/1674-4527/15/8/002
  * Maíz Apellániz & Weiler (2018) Maíz Apellániz, J., & Weiler, M. 2018, A&A, 619, A180, doi: 10.1051/0004-6361/201834051
  * Onken et al. (2019) Onken, C. A., Wolf, C., Bessell, M. S., et al. 2019, PASA, 36, e033, doi: 10.1017/pasa.2019.27
  * Peek & Graves (2010) Peek, J. E. G., & Graves, G. J. 2010, ApJ, 719, 415, doi: 10.1088/0004-637X/719/1/415
  * Riello et al. (2018) Riello, M., De Angeli, F., Evans, D. W., et al. 2018, A&A, 616, A3, doi: 10.1051/0004-6361/201832712
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Weiler (2018) Weiler, M. 2018, A&A, 617, A138, doi: 10.1051/0004-6361/201833462
  * Wu et al. (2011) Wu, Y., Luo, A. L., Li, H.-N., et al. 2011, Research in Astronomy and Astrophysics, 11, 924, doi: 10.1088/1674-4527/11/8/006
  * Wu et al. (2019) Wu, Y., Xiang, M., Zhao, G., et al. 2019, MNRAS, 484, 5315, doi: 10.1093/mnras/stz256
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Yuan et al. (2015a) Yuan, H., Liu, X., Xiang, M., Huang, Y., & Chen, B. 2015a, ApJ, 799, 134, doi: 10.1088/0004-637X/799/2/134
  * Yuan et al. (2015b) Yuan, H., Liu, X., Xiang, M., et al. 2015b, ApJ, 799, 135, doi: 10.1088/0004-637X/799/2/135
  * Yuan et al. (2015c) —. 2015c, ApJ, 799, 133, doi: 10.1088/0004-637X/799/2/133
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723, doi: 10.1088/1674-4527/12/7/002


Triangle Groups: Automorphic Forms and Nonlinear Differential Equations

Triangle Groups: Automorphic Forms
and Nonlinear Differential Equations

Sujay K. ASHOK †, Dileep P. JATKAR ‡ and Madhusudhan RAMAN §

S.K. Ashok, D.P. Jatkar and M. Raman

† Institute of Mathematical Sciences, Homi Bhabha National Institute (HBNI),
† IV Cross Road, C.I.T. Campus, Taramani, Chennai 600 113, India
<EMAIL_ADDRESS>
‡ Harish-Chandra Research Institute, Homi Bhabha National Institute (HBNI),
‡ Chhatnag Road, Jhunsi, Allahabad 211 019, India<EMAIL_ADDRESS>
§ Department of Theoretical Physics, Tata Institute of Fundamental Research,
§ Homi Bhabha Road, Navy Nagar, Colaba, Mumbai 400 005, India
<EMAIL_ADDRESS>
Received April 21, 2020, in final form October 05, 2020; Published online
October 11, 2020

We study the relations governing the ring of quasiautomorphic forms associated
to triangle groups with a single cusp, thereby extending our earlier results
on Hecke groups. The Eisenstein series associated to these triangle groups are
shown to satisfy Ramanujan-like identities. These identities in turn allow us
to associate a nonlinear differential equation to each triangle group. We show
that they are solved by the quasiautomorphic weight-2 Eisenstein series
associated to the triangle group and its orbit under the group action. We
conclude by discussing the Painlevé property of these nonlinear differential
equations.

triangle groups; Chazy equations; Painlevé analysis

34M55; 11F12; 33E30

##  1 Introduction and discussion

Triangle groups are an infinite family of discrete subgroups of ${\rm
PSL}(2,\mathbb{R})$, the group of orientation-preserving isometries of the
upper-half plane. Abstractly, these groups are characterised by an ordered
triple of positive integers $\\{m_{i}\\}_{i=1}^{3}$ and generated by the
letters $\\{g_{i}\\}_{i=1}^{3}$, and circumscribed by the relations

$\displaystyle g_{i}^{m_{i}}=1\qquad\text{for}\
i\in\\{1,2,3\\}\qquad\text{and}\qquad g_{1}g_{2}g_{3}=1.$

The theory of automorphic functions and forms associated to triangle groups
has been detailed in [11]. Triangle groups with arithmetic properties have
appeared in the study of string theories [7, 14, 19], supersymmetric gauge
theories [3], the study of knots [17, 25, 26], and simple quantum mechanical
systems [5, 23] and further motivates our investigations. The present work can
be thought of as an extension of our results for Hecke groups obtained in [4]
to all triangle group with a single cusp.

Differential equations satisfied by modular [13, 16] and quasi-modular [1, 6,
21] objects have been the subject of sustained interest. The primary objects
of study in this paper are an infinite class of nonlinear differential
equations naturally associated to triangle groups. In Section 2, we review the
construction of quasi-automorphic forms associated to arbitrary triangle
groups as solutions to the generalised Halphen system, a system of first-order
differential equations, following [11]. The automorphic forms associated to
these triangle groups are then found to satisfy a system of first-order
differential equations whose structure is strongly reminiscent of the
Ramanujan identities governing the ring of (quasi-)modular forms of ${\rm
PSL}(2,\mathbb{Z})$.

The construction of automorphic forms also exhibits the existence of certain
ring relations – that is, a set of relations that imply that the ring of
automorphic forms associated to a triangle group is not freely generated.
These ring relations play a central role in this note. Indeed, in Section 3 we
use both the Ramanujan identities and the ring relations to derive a non-
linear differential equation previously derived in [15]. Notably, our
derivation does not rely on any explicit Fourier expansions or require the use
of computer algebra.

In the final section, we present the derivation of higher-order differential
equations associated to each triangle group via differential elimination on
the Ramanujan identities. These equations are natural analogues of the Chazy
equation associated to ${\rm PSL}(2,\mathbb{Z})$. We then prove that these
differential equations possess the Painlevé property, which is to say that the
only movable singularities of these differential equations are poles.

In his analysis of second-order ordinary differential equations, Painlevé
restricted his attention to only those equations whose movable singularities
were poles. All these equations turned out to be integrable. Based on this,
the authors of [2] devised the Painlevé test, a set of criteria that, if
satisfied by a differential equation, establishes the Painlevé property. This
test works perfectly for all soliton equations but requires a more careful
stability analysis for those equations which possess negative resonances,
along the lines of [9]. The nonlinear differential equations we study in this
note generically possess negative resonances. In some cases, as we will show,
this is no obstruction to the demonstration of the Painlevé property.

A proof of the Painlevé property is a prerequisite to the demonstration of
integrability. Indeed, the higher-order Chazy equations uncovered in [4]
already possess tantalising hints that these equations may be integrable.
Consider as an example the Chazy equation corresponding to the Hecke group
${\rm H}(4)$ in [4, Theorem 4.3], written in terms of the Ramanujan–Serre
derivative D and a weight-$4$ modular form $Z$:

$\displaystyle\mathrm{D}^{3}Z+6Z\mathrm{D}Z=0.$

This presentation manifests a striking resemblance to the stationary
Korteweg–de Vries equation:

$\displaystyle u^{\prime\prime\prime}+6uu^{\prime}=0.$

Given the uniform structure of the Ramanujan identities for triangle groups we
write down, and the higher-order Chazy and Maier equations that are consequent
upon them, one is tempted to look for a framework, perhaps even a hierarchy,
within which these differential equations may be organised. The results of
this note, in particular investigations into the Painlevé property, are an
important step in this direction.

##  2 Triangle groups and Ramanujan identities

We begin with a brief review of triangle groups and the associated generalised
Halphen system. We will then review the main theorems of [11] regarding the
construction of Eisenstein series in terms of the solutions to the Halphen
system. Starting from these definitions, we then derive the Ramanujan
relations that are valid for any triangle group with a single cusp.

###  2.1 Triangle groups

A triangle group is identified by its _type_ , denoted by $\mathfrak{t}$ and
consisting of a triple of integers $2\leq m_{1}\leq m_{2}\leq
m_{3}\leq\infty$. The integers that make up a triple define a triple of angles
$(\pi/m_{1},\pi/m_{2},\pi/m_{3})$, which are taken to be the angles subtended
by the edges of a triangle in the hyperbolic plane. The group of reflections
of this triangle across its sides generates a tessellation of the hyperbolic
plane, and is called a triangle group.

Triangle groups are discrete subgroups of ${\rm PSL}(2,\mathbb{R})$, the group
of isometries of the upper-half plane $\mathbb{H}$, and we study those whose
fundamental domains have finite hyperbolic area.111As a result, the following
types are not considered: (i) $(2,2,m)$ for all $m\leq\infty$, (ii) $(2,3,n)$
for $n\leq 6$, and (iii) the types $(2,4,4)$ and $(3,3,3)$. We will further
restrict our attention to those triangle groups with a single cusp – meaning
the angle subtended by one pair of edges of the aforementioned triangle is
zero – and without loss of generality this choice will be implemented by
setting $m_{3}=\infty$.222For a more algebraic characterisation, these
triangle groups are isomorphic to the free product group
$\mathbb{Z}_{m_{1}}*\mathbb{Z}_{m_{2}}$, where the notation $\mathbb{Z}_{k}$
denotes a cyclic group of order $k$. In what follows, when we refer to
triangle groups we will always mean triangle groups with a single cusp,
characterised by a pair of integers so its type
$\mathfrak{t}=(m_{1},m_{2},\infty)$ such that $2\leq m_{1}\leq m_{2}<\infty$.
For example, the modular group ${\rm PSL}(2,\mathbb{Z})$ is a triangle group
of type $(2,3,\infty)$, and the Hecke groups ${\rm H}(m)$ are of type
$(2,m,\infty)$. Finally, throughout this paper, the subscript $\mathfrak{t}$
will indicate that the object in question is associated to any specific choice
of triangle group.

We now move on to discuss the construction of Eisenstein series associated to
these triangle groups. This will be done following [11], who construct these
automorphic forms out of solutions to a generalised Halphen system.

###  2.2 Generalized Halphen systems and Eisenstein series

The generalized Halphen system is a set of coupled ordinary first-order
differential equations for three variables
$\\{t_{k,\mathfrak{t}}(\tau)\\}_{k=1}^{3}$ given as follows:

$\displaystyle
t_{1,\mathfrak{t}}^{\prime}=(a-1)(t_{1,\mathfrak{t}}t_{2,\mathfrak{t}}+t_{1,\mathfrak{t}}t_{3,\mathfrak{t}}-t_{2,\mathfrak{t}}t_{3,\mathfrak{t}})+(b+c-1)t_{1,\mathfrak{t}}^{2},$
$\displaystyle
t_{2,\mathfrak{t}}^{\prime}=(b-1)(t_{2,\mathfrak{t}}t_{3,\mathfrak{t}}+t_{2,\mathfrak{t}}t_{1,\mathfrak{t}}-t_{1,\mathfrak{t}}t_{3,\mathfrak{t}})+(a+c-1)t_{2,\mathfrak{t}}^{2},$
$\displaystyle
t_{3,\mathfrak{t}}^{\prime}=(c-1)(t_{3,\mathfrak{t}}t_{1,\mathfrak{t}}+t_{3,\mathfrak{t}}t_{2,\mathfrak{t}}-t_{1,\mathfrak{t}}t_{2,\mathfrak{t}})+(a+b-1)t_{3,\mathfrak{t}}^{2}.$
(2.1)

Here, the parameters $\\{a,b,c\\}$ are specified by the type $\mathfrak{t}$ of
the triangle group as

$\displaystyle
a=\frac{1}{2}\left(1-\frac{1}{m_{1}}+\frac{1}{m_{2}}\right)\\!,\qquad
b=\frac{1}{2}\left(1-\frac{1}{m_{1}}-\frac{1}{m_{2}}\right)\\!,\qquad
c=\frac{1}{2}\left(1+\frac{1}{m_{1}}-\frac{1}{m_{2}}\right)\\!,$ (2.2)

and the accent ′ denotes the following derivative:

$\displaystyle
h^{\prime}:=\frac{1}{2\pi\text{i}}\frac{\text{d}}{\text{d}\tau}h(\tau),$

where $\tau$ is a coordinate on the upper-half plane $\mathbb{H}$.

The solutions to the generalized Halphen system can be obtained explicitly in
terms of hypergeometric functions [11, Theorem 3(i)]. We relegate a discussion
of these solutions to Appendix A – notably, our conventions differ slightly
from theirs – and instead proceed to the construction of Eisenstein series. We
start with the following theorem due to [11, see p. 707 and Theorem 4(iv)]:

######  Theorem 2.1.

The algebra of holomorphic automorphic forms of a triangle group
$\mathfrak{t}$ is generated by two branches of Eisenstein series:
$\big{\\{}E_{2k,\mathfrak{t}}^{(1)}\big{\\}}_{k=2}^{m_{1}}$ and
$\big{\\{}E_{2k,\mathfrak{t}}^{(2)}\big{\\}}_{k=2}^{m_{2}}$. For each
admissible $k$, the corresponding Eisenstein series has an automorphic weight
$2k$ under the triangle group. These Eisenstein series are constructed out of
solutions to the generalised Halphen system as follows: define the linear
combinations

$\displaystyle\mathsf{x_{\mathfrak{t}}}=t_{1,\mathfrak{t}}-t_{2,\mathfrak{t}}\qquad\text{and}\qquad\mathsf{y_{\mathfrak{t}}}=t_{3,\mathfrak{t}}-t_{2,\mathfrak{t}},$
(2.3)

in terms of which the Eisenstein series are

$\displaystyle
E_{2k,\mathfrak{t}}^{(1)}=\mathsf{x}_{\mathfrak{t}}\mathsf{y}_{\mathfrak{t}}^{k-1}\qquad$
$\displaystyle\text{for}\ 2\leq k\leq m_{1},$ (2.4) $\displaystyle
E_{2k,\mathfrak{t}}^{(2)}=\mathsf{x}_{\mathfrak{t}}^{k-1}\mathsf{y}_{\mathfrak{t}}\qquad$
$\displaystyle\text{for}\ 2\leq k\leq m_{2}.$ (2.5)

The above construction of the Eisenstein series in terms of solutions to the
generalised Halphen system straightforwardly implies the following results.

######  Corollary 2.2.

The algebra of automorphic forms is not freely generated. For both branches of
Eisenstein series, i.e., for $\ell\in\\{1,2\\}$ we have the following
relations:

$E_{2p,\mathfrak{t}}^{(\ell)}E_{2(k-p+1),\mathfrak{t}}^{(\ell)}=E_{2p^{\prime},\mathfrak{t}}^{(\ell)}E_{2(k-p^{\prime}+1),\mathfrak{t}}^{(\ell)},$
(2.6)

for any integers $2\leq p\leq k-1$ and $2\leq p^{\prime}\leq k-1$.
Additionally, we have the following relations binding the two branches of
Eisenstein series:

$E_{2k,\mathfrak{t}}^{(1)}E_{2k,\mathfrak{t}}^{(2)}=(E_{4,\mathfrak{t}})^{k}.$
(2.7)

Collectively, we will refer to the above identities as ring relations. They
will be used extensively in the following sections.

######  Remark 2.3.

Associated to every triangle group, there exists a unique weight-$4$
Eisenstein series, obtained as the case $k=2$ of (2.5). Since this automorphic
form may be thought of as belonging to either branch of Eisenstein series, we
denote it without any reference to these branches, as $E_{4,\mathfrak{t}}$.

The triangle groups also come equipped with a quasi-automorphic weight $2$
Eisenstein series $E_{2,\mathfrak{t}}$, which is also given in terms of a
linear combination of the generalised Halphen variables $t_{i,\mathfrak{t}}$
[11, Theorem 4(iii)].

######  Theorem 2.4.

For any triangle group with type $\mathfrak{t}$, we can associate a
holomorphic, quasi-automorphic weight-$2$ Eisenstein series
$E_{2,\mathfrak{t}}$. In terms of solutions to a generalised Halphen system,
it is given by

$\displaystyle
E_{2,\mathfrak{t}}=\frac{1}{m_{1}+m_{2}-m_{1}m_{2}}\big{(}2m_{1}\mathsf{x}_{\mathfrak{t}}+2m_{2}\mathsf{y}_{\mathfrak{t}}+(m_{1}+m_{2}+m_{1}m_{2})t_{2,\mathfrak{t}}\big{)}.$
(2.8)

This choice of linear combination is tied to the choice of normalisation for
solutions to the generalised Halphen system and ensures that at
$\text{i}\infty$, $E_{2,\mathfrak{t}}$ is unity – see Appendix A.

With the introduction of the Eisenstein series that generate the ring of
holomorphic, quasi-automophic forms of triangle groups $\mathfrak{t}$ now
complete, we are in a position to state our first result. It generalizes our
previous results on Ramanujan identities for Hecke groups [4, 22] to an
arbitrary triangle group with cusp.

######  Lemma 2.5.

The Eisenstein series associated to a triangle group $\mathfrak{t}$, as
defined in Theorems 2.1 and 2.4, and for $k>2$, obey the following identities:

$\displaystyle
E_{2,\mathfrak{t}}^{\prime}=\frac{(m_{1}m_{2}-m_{1}-m_{2})}{2m_{1}m_{2}}\left(E_{2,\mathfrak{t}}^{2}-E_{4,\mathfrak{t}}\right)\\!,$
$\displaystyle
E_{4,\mathfrak{t}}^{\prime}=2\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}E_{4,\mathfrak{t}}-\frac{m_{1}-2}{m_{1}}E_{6,\mathfrak{t}}^{(1)}-\frac{m_{2}-2}{m_{2}}E_{6,\mathfrak{t}}^{(2)},$
$\displaystyle
E_{2k,\mathfrak{t}}^{{(1)}^{\prime}}=\frac{k(m_{1}m_{2}\\!-\\!m_{1}\\!-\\!m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}E_{2k,\mathfrak{t}}^{(1)}\\!-\\!\frac{(km_{2}\\!-\\!k\\!-\\!m_{2})}{m_{2}}E_{4,\mathfrak{t}}E_{2k\\!-\\!2,\mathfrak{t}}^{(1)}\\!-\\!\frac{(m_{1}\\!-\\!k)}{m_{1}}E_{2k+2,\mathfrak{t}}^{(1)},$
$\displaystyle
E_{2k,\mathfrak{t}}^{{(2)}^{\prime}}=\frac{k(m_{1}m_{2}\\!-\\!m_{1}\\!-\\!m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}E_{2k,\mathfrak{t}}^{(2)}\\!-\\!\frac{(km_{1}\\!-\\!k\\!-\\!m_{1})}{m_{1}}E_{4,\mathfrak{t}}E_{2k\\!-\\!2,\mathfrak{t}}^{(2)}\\!-\\!\frac{(m_{2}\\!-\\!k)}{m_{2}}E_{2k+2,\mathfrak{t}}^{(2)}.$
(2.9)

###### Proof.

We perform a linear transformation from the generalized Halphen variables
$\\{t_{i,\mathfrak{t}}\\}_{i=1}^{3}$ to the variables
$(E_{2,\mathfrak{t}},\mathsf{x}_{\mathfrak{t}},\mathsf{y}_{\mathfrak{t}})$
using (2.3) and (2.8). The generalised Halphen system is then equivalent to
the following system of ordinary first-order differential equations:

$\displaystyle
E_{2,\mathfrak{t}}^{\prime}=\frac{(m_{1}m_{2}-m_{1}-m_{2})}{2m_{1}m_{2}}\left(E_{2,\mathfrak{t}}^{2}-E_{4,\mathfrak{t}}\right)\\!,$
$\displaystyle\mathsf{x}_{\mathfrak{t}}^{\prime}=\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}\mathsf{x}_{\mathfrak{t}}+\frac{\mathsf{x}_{\mathfrak{t}}^{2}}{m_{2}}-\left(1-\frac{1}{m_{1}}\right)E_{4,\mathfrak{t}},$
$\displaystyle\mathsf{y}_{\mathfrak{t}}^{\prime}=\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}\mathsf{y}_{\mathfrak{t}}+\frac{\mathsf{y}_{\mathfrak{t}}^{2}}{m_{1}}-\left(1-\frac{1}{m_{2}}\right)E_{4,\mathfrak{t}}.$
(2.10)

The first equation in (2.10) proves the first of the Ramanujan identities. For
$E_{4,\mathfrak{t}}$, we consider the derivative of (2.5) for $k=2$ to get

$\displaystyle
E_{4,\mathfrak{t}}^{\prime}=\left[\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}\mathsf{x}_{\mathfrak{t}}+\frac{\mathsf{x}_{\mathfrak{t}}^{2}}{m_{2}}-\left(1-\frac{1}{m_{1}}\right)\mathsf{x}_{\mathfrak{t}}\mathsf{y}_{\mathfrak{t}}\right]\mathsf{y}_{\mathfrak{t}}$
$\displaystyle\hphantom{E_{4,\mathfrak{t}}^{\prime}=}{}+\mathsf{x}_{\mathfrak{t}}\left[\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}\mathsf{y}_{\mathfrak{t}}+\frac{\mathsf{y}_{\mathfrak{t}}^{2}}{m_{1}}-\left(1-\frac{1}{m_{2}}\right)\mathsf{x}_{\mathfrak{t}}\mathsf{y}_{\mathfrak{t}}\right]$
$\displaystyle\hphantom{E_{4,\mathfrak{t}}^{\prime}}{}=2\frac{(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}m_{2}}E_{2,\mathfrak{t}}E_{4,\mathfrak{t}}-\frac{m_{1}-2}{m_{1}}E_{6,\mathfrak{t}}^{(1)}-\frac{m_{2}-2}{m_{2}}E_{6,\mathfrak{t}}^{(2)}.$

Similarly, for the relations governing either branch of Eisenstein series
$E_{2k,\mathfrak{t}}^{(\ell)}$ of weight $2k>4$, differentiate (2.5) and use
(2.10), then resolve all products of $\mathsf{x_{\mathfrak{t}}}$ and
$\mathsf{y_{\mathfrak{t}}}$ into Eisenstein series according to (2.5). ∎

Owing to their similarities with relations satisfied by the Eisenstein series
associated to $\mathrm{PSL}(2,\mathbb{Z})$, these identities will henceforth
be referred to as the Ramanujan identities associated to the triangle group
$\mathfrak{t}$.

##  3 Maier equations for triangle groups

The Ramanujan identities are a set of coupled first order differential
equations that involve all the Eisenstein series. However as we have seen,
there are non-trivial ring relations that relate the higher weight Eisenstein
series. We now show how it is possible to use the ring relations and deduce
from the Ramanujan relations a higher order non-linear differential equation
for every triangle group. We show that the equation coincides precisely with
the one obtained by Maier [15]. Our derivation has the advantage of being
completely elementary without the need for advanced computer algebra or
explicit Fourier expansions of the automorphic forms and turns out to follow
simply from the existence of the Ramanujan identities and the ring relations.

Before we proceed further we first introduce the notion of a covariant
derivative. In fact, the structure of the Ramanujan identities lends itself to
a natural choice of derivation on ring of automorphic forms. Let the space of
weight-$k$ automorphic forms associated to a triangle group $\mathfrak{t}$ be
denoted $\mathfrak{m}_{k}$.333The dependence on the choice of triangle group
is left implicit to avoid clutter. We define the operator D such that

$\displaystyle\text{D}\colon\
\mathfrak{m}_{k}\longrightarrow\mathfrak{m}_{k+2},$

Explicitly, this Ramanujan–Serre derivative is defined as

$\displaystyle\text{D}=\frac{1}{2\pi\text{i}}\frac{\text{d}missing}{\text{d}\tau}-\frac{k}{2}\left(1-\frac{1}{m_{1}}-\frac{1}{m_{2}}\right)E_{2,\mathfrak{t}}.$

With this definition, the first few Ramanujan identities may be written in
terms of the Ramanujan–Serre derivative as

$\displaystyle\text{D}E_{4,\mathfrak{t}}=\frac{2-m_{1}}{m_{1}}E_{6,\mathfrak{t}}^{(1)}+\frac{2-m_{2}}{m_{2}}E_{6,\mathfrak{t}}^{(2)},$
(3.1)
$\displaystyle\text{D}E_{6,\mathfrak{t}}^{(1)}=\frac{3-m_{1}}{m_{1}}E_{8,\mathfrak{t}}^{(1)}+\frac{3-2m_{2}}{m_{2}}E_{4,\mathfrak{t}}^{2},$
(3.2)
$\displaystyle\text{D}E_{6,\mathfrak{t}}^{(2)}=\frac{3-m_{2}}{m_{2}}E_{8,\mathfrak{t}}^{(2)}+\frac{3-2m_{1}}{m_{1}}E_{4,\mathfrak{t}}^{2}.$
(3.3)

We now act with the Serre derivative on (3.1), and use equations (3.2) and
(3.3) along with the ring relation (2.6), i.e., we use

$\displaystyle
E_{4,\mathfrak{t}}^{(\ell)}E_{8,\mathfrak{t}}^{(\ell)}=\big{(}E_{6,\mathfrak{t}}^{(\ell)}\big{)}^{2}\qquad\text{for}\
\ell=1,2,$

to obtain the following differential equation relating Eisenstein series’ of
weights four and six:

$\displaystyle
E_{4,\mathfrak{t}}\text{D}^{2}E_{4,\mathfrak{t}}=\frac{(3-2m_{1})(2-m_{2})+(3-2m_{2})(2-m_{1})}{m_{1}m_{2}}E_{4,\mathfrak{t}}^{3}$
$\displaystyle\hphantom{E_{4,\mathfrak{t}}\text{D}^{2}E_{4,\mathfrak{t}}=}{}+\frac{(m_{1}-2)(m_{1}-3)}{m_{1}^{2}}\big{(}E_{6,\mathfrak{t}}^{(1)}\big{)}^{2}+\frac{(m_{2}-2)(m_{2}-3)}{m_{2}^{2}}\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}.$
(3.4)

Note that each term in the above equation is automorphic with weight $12$, and
consequently that the equation is invariant under triangle group action.

Before turning to the main goal of this section, i.e., deriving the Maier
equations for arbitrary triangle groups [15], we first consider two special
cases. These cases are distinguished by simpler Maier equations with lower
weight than one might expect, chiefly due to the simplifications that arise
when the integers specifying the type $\mathfrak{t}$ of the triangle group are
tuned to special values.

######  Theorem 3.1 (isosceles triangle groups).

The Maier equations corresponding to isosceles triangle groups with type
$\mathfrak{t}_{M}=(M,M,\infty)$ are identical to the Maier equations for the
Hecke group $\mathrm{H}(M)$.

###### Proof.

Equations (3.1) and (3.4) simplify considerably for isosceles triangle groups
with $m_{1}=m_{2}=M$. Using the square of (3.1) and the ring relation (2.7),
i.e.,

$\displaystyle
E_{6,\mathfrak{t}_{M}}^{(1)}E_{6,\mathfrak{t}_{M}}^{(2)}=E_{4,\mathfrak{t}_{M}}^{3},$
(3.5)

we can write (3.4) as

$\displaystyle
E_{4,\mathfrak{t}_{M}}\text{D}^{2}E_{4,\mathfrak{t}_{M}}=\frac{M-3}{M-2}(\text{D}E_{4,\mathfrak{t}_{M}})^{2}+\frac{2(M-2)}{M}E_{4,\mathfrak{t}_{M}}^{3}.$

Up to an unimportant rescaling of $E_{4,\mathfrak{t}_{M}}$, this is the same
equation that was derived for the Hecke group ${\rm H}(M)$ in [4]. ∎

This fact was already observed in [15] and we see that it is a simple
consequence of the Ramanujan relations and the ring relations. There exists a
simple geometrical basis for this fact: an isosceles hyperbolic triangle of
type $(M,M,\infty)$ can be bisected into two hyperbolic triangles of type
$(2,M,\infty)$ [18]. We now turn to triangles groups of type
$\mathfrak{t}_{3,M}=(3,M,\infty)$. We present this example because it is
instructive and will help build intuition for the most general case.

######  Example 3.2 (type $(3,M,\infty)$ triangle groups).

In this case as well the goal is to write an equation purely in terms of the
modular covariant derivatives $\text{D}^{k}E_{4,\mathfrak{t}_{3,M}}$. Start
with (3.1) and solve for $E_{6,\mathfrak{t}_{3,M}}^{(1)}$ as

$\displaystyle
E_{6,\mathfrak{t}_{3,M}}^{(1)}=-3\text{D}E_{4,\mathfrak{t}_{3,M}}-\frac{3(M-2)}{M}E_{6,\mathfrak{t}_{3,M}}^{(2)}.$

By combining this with the ring relation (3.5), we obtain the relation:

$\displaystyle\text{D}E_{4,\mathfrak{t}_{3,M}}E_{6,\mathfrak{t}_{3,M}}^{(2)}=-\frac{M-2}{M}\big{(}E_{6,\mathfrak{t}_{3,M}}^{(2)}\big{)}^{2}-\frac{1}{3}E_{4,\mathfrak{t}_{3,M}}^{3}.$

One can write an independent expression for the square of
$E_{6,\mathfrak{t}_{3,M}}^{(2)}$ in terms of the $E_{4,\mathfrak{t}_{3,M}}$
and its modular covariant derivatives courtesy of (3.4):

$\displaystyle\big{(}E_{6,\mathfrak{t}_{3,M}}^{(2)}\big{)}^{2}=\frac{M^{2}}{(M-3)(M-2)}E_{4,\mathfrak{t}_{3,M}}\text{D}^{2}E_{4,\mathfrak{t}_{3,M}}+\frac{M(5M-9)}{3(M-3)(M-2)}(E_{4,\mathfrak{t}_{3,M}})^{3}.$
(3.6)

This leads to the following relation:

$\displaystyle\text{D}E_{4,\mathfrak{t}_{3,M}}E_{6,\mathfrak{t}_{3,M}}^{(2)}=\frac{2}{3}\frac{(2M-3)}{(M-3)}(E_{4,\mathfrak{t}_{3,M}})^{3}-\frac{M}{M-3}E_{4,\mathfrak{t}_{3,M}}\text{D}^{2}E_{4,\mathfrak{t}_{3,M}}.$

Squaring both sides of this equation and using the ring relation (3.6) once
again leads to the following equation, with each term of weight 20:

$\displaystyle
9M^{2}(M-3)\text{D}^{2}E_{4,\mathfrak{t}_{3,M}}(\text{D}E_{4,\mathfrak{t}_{3,M}})^{2}-9M^{2}(M-2)(\text{D}^{2}E_{4,\mathfrak{t}_{3,M}})^{2}E_{4,\mathfrak{t}_{3,M}}$
$\displaystyle\qquad{}+3M(M-3)(9-5M)(\text{D}E_{4,\mathfrak{t}_{3,M}})^{2}(E_{4,\mathfrak{t}_{3,M}})^{2}$
$\displaystyle\qquad{}+12M(M-2)(2M-3)\text{D}^{2}E_{4,\mathfrak{t}_{3,M}}(E_{4,\mathfrak{t}_{3,M}})^{3}-4(M-2)(2M-3)^{2}(E_{4,\mathfrak{t}_{3,M}})^{5}=0.$

We now bring this strategy to bear on arbitrary triangle groups. This theorem
corresponds directly to [15, Theorem 6.4] and we use the same notation for the
coefficients for ease of comparison.

######  Theorem 3.3.

For an arbitrary triangle group of type $\mathfrak{t}=(m_{1},m_{2},\infty)$,
the weight four automorphic form $E_{4,\mathfrak{t}}$ satisfies the following
weight-$24$ differential equation:

$\displaystyle
C_{88}E_{4,\mathfrak{t}}^{2}\big{(}\mathrm{D}^{2}E_{4,\mathfrak{t}}\big{)}^{2}+C_{86}E_{4,\mathfrak{t}}\big{(}\mathrm{D}E_{4,\mathfrak{t}}\big{)}^{2}\mathrm{D}^{2}E_{4,\mathfrak{t}}+C_{84}E_{4,\mathfrak{t}}^{4}\big{(}\mathrm{D}^{2}E_{4,\mathfrak{t}}\big{)}+C_{66}(\mathrm{D}E_{4,\mathfrak{t}})^{4}$
$\displaystyle\hphantom{C_{88}E_{4,\mathfrak{t}}^{2}\big{(}\mathrm{D}^{2}E_{4,\mathfrak{t}}\big{)}^{2}}{}+C_{64}E_{4,\mathfrak{t}}^{3}(\mathrm{D}E_{4,\mathfrak{t}})^{2}+C_{44}E_{4,\mathfrak{t}}^{6}=0,$
(3.7)

where the coefficients are given by

$\displaystyle C_{88}=m_{1}^{2}m_{2}^{2}(m_{1}-2)(m_{2}-2),$ $\displaystyle
C_{86}=-m_{1}^{2}m_{2}^{2}((m_{1}-2)(m_{2}-3)+(m_{2}-2)(m_{1}-3)),$
$\displaystyle C_{84}=-4m_{1}m_{2}(m_{1}-2)(m_{2}-2)(m_{1}m_{2}-m_{1}-m_{2}),$
$\displaystyle C_{66}=m_{1}^{2}m_{2}^{2}(m_{1}-3)(m_{2}-3),$ $\displaystyle
C_{64}=m_{1}m_{2}\left(3(m_{1}-m_{2})^{2}+2m_{2}(m_{2}-3)(m_{1}-2)^{2}+2m_{1}(m_{1}-3)(m_{2}-2)^{2}\right)\\!,$
$\displaystyle C_{44}=4(m_{1}-2)(m_{2}-2)(m_{1}m_{2}-m_{1}-m_{2})^{2}.$ (3.8)

###### Proof.

First, solve for $E_{6,\mathfrak{t}}^{(1)}$ using (3.1):

$\displaystyle
E_{6,\mathfrak{t}}^{(1)}=-\frac{\text{D}E_{4,\mathfrak{t}}m_{1}}{m_{1}-2}-\frac{m_{1}\left(m_{2}-2\right)}{m_{2}\left(m_{1}-2\right)}E_{6,\mathfrak{t}}^{(2)}.$

Then substitute this into the ring relation (3.5) to obtain the following
equation:

$\displaystyle\frac{m_{2}(m_{1}-2)}{m_{1}(m_{2}-2)}E_{4,\mathfrak{t}}^{3}+\frac{m_{2}}{m_{2}-2}E_{6,\mathfrak{t}}^{(2)}\text{D}E_{4,\mathfrak{t}}+\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}=0.$
(3.9)

A second independent weight-12 equation is arrived at by acting with
$E_{4,\mathfrak{t}}\text{D}$ on (3.1) and in the resulting equation, we
substitute the equations (3.2) and (3.3). After this, on taking suitable
linear combinations we eliminate $\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}$
and obtain the following equation:

$\displaystyle\text{D}E_{4,\mathfrak{t}}E_{6,\mathfrak{t}}^{(2)}=-\frac{(m_{1}-2)m_{2}}{m_{1}-m_{2}}E_{4,\mathfrak{t}}\text{D}^{2}E_{4,\mathfrak{t}}+\frac{(m_{1}-3)m_{2}}{m_{1}-m_{2}}(\text{D}E_{4,\mathfrak{t}})^{2}$
$\displaystyle\hphantom{\text{D}E_{4,\mathfrak{t}}E_{6,\mathfrak{t}}^{(2)}=}{}+\frac{2(m_{1}-2)(m_{1}m_{2}-m_{1}-m_{2})}{m_{1}(m_{1}-m_{2})}E_{4,\mathfrak{t}}^{3}.$
(3.10)

We substitute this into (3.9) to obtain the ring relation for
$\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}$:

$\displaystyle\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}=-\frac{(m_{1}-2)m_{2}^{2}}{(m_{1}-m_{2})(m_{2}-2)}E_{4,\mathfrak{t}}\text{D}^{2}E_{4,\mathfrak{t}}+\frac{(m_{1}-3)m_{2}^{2}}{(m_{1}-m_{2})(m_{2}-2)}\text{D}E_{4,\mathfrak{t}}^{2}$
$\displaystyle\hphantom{\big{(}E_{6,\mathfrak{t}}^{(2)}\big{)}^{2}=}{}+\frac{m_{2}(m_{1}-2)(2m_{1}m_{2}-3m_{1}-m_{2})}{m_{1}(m_{1}-m_{2})(m_{2}-2)}E_{4,\mathfrak{t}}^{3}.$
(3.11)

We now have all the ingredients in place to derive the Maier equation for the
general triangle group. Consider the square of (3.10) and substitute the
expression for the square of $E_{6,\mathfrak{t}}^{(2)}$ derived in (3.11)
above. This leads to the weight-$24$ equation written purely in terms of
$E_{4,\mathfrak{t}}$ and its modular covariant derivatives, given in (3.7),
with the coefficients given in (3.8). Finally, that the orbit of
$E_{4,\mathfrak{t}}$ under the corresponding triangle group also solves the
nonlinear differential equation under consideration follows from its manifest
covariance. ∎

Up to a normalization factor, this precisely matches the Maier equation for
the triangle group of type $(m_{1},m_{2},\infty)$ that is presented in [15,
Theorem 6.4]. We see that this differential equation is a direct consequence
of the Ramanujan identities satisfied by the Eisenstein series in combination
with the algebraic ring relations that bind them.

##  4 Painlevé analysis

We now turn to a study of the nonlinear differential equations that we
introduced in the previous section. In particular, we will show that these
equations satisfy the Painlevé property, i.e., that the only movable
singularities of these differential equations are poles. For a brief review of
the stability of differential equations and the Painlevé property, we refer
the reader to [4].

###  4.1 Chazy equations

In this section, we outline a process of differential elimination that may be
used to derive higher-order nonlinear differential equations satisfied by
$E_{2,\mathfrak{t}}$. These differential equations are of order $m_{2}$, and
are natural analogues of the higher-order Chazy equations studied in [4].444We
remind the reader that without loss of generality, we choose the integers that
specify the type $\mathfrak{t}$ of the triangle group such that $m_{1}\leq
m_{2}$.

The process of differential elimination is complicated by presence of two
distinct branches of Eisenstein series. In order to sidestep this difficulty,
we effect a change of basis in the space of automorphic forms that is
“diagonal” in the sense that only specific linear combinations of the two
branches of Eisenstein series appear together in the Ramanujan identities. For
triangle group $\mathfrak{t}$ and $k\geq 3$, define the following combinations
of the two branches of Eisenstein series:

$\displaystyle
G_{2k,\mathfrak{t}}=\prod_{p=2}^{k-1}\left(\frac{m_{1}-p}{m_{1}}\right)E_{2k,\mathfrak{t}}^{(1)}+\prod_{p=2}^{k-1}\left(\frac{m_{2}-p}{m_{2}}\right)E_{2k,\mathfrak{t}}^{(2)}.$
(4.1)

Since the Eisenstein series’ of weights two and four are unique we identify
$G_{4,\mathfrak{t}}=E_{4,\mathfrak{t}}$ and
$G_{2,\mathfrak{t}}=E_{2,\mathfrak{t}}$ for uniformity of notation. Once this
change of basis is effected, we may proceed with straightforward differential
elimination. Note that the $G_{2k,\mathfrak{t}}$ are not algebraically
independent. This is once again a consequence of the ring relations. The
utility of this change of basis is presented in the following simple example.

######  Example 4.1 (type $(4,4,\infty)$ triangle group).

In terms of the linear combinations defined in (4.1), the Ramanujan identities
(2.9) associated to the triangle group $(4,4,\infty)$ take the following form:

$\displaystyle
G_{2,\mathfrak{t}}^{\prime}=\frac{1}{4}\big{(}G_{2,\mathfrak{t}}^{2}-G_{4,\mathfrak{t}}\big{)},$
$\displaystyle
G_{4,\mathfrak{t}}^{\prime}=G_{2,\mathfrak{t}}G_{4,\mathfrak{t}}-G_{6,\mathfrak{t}},$
$\displaystyle
G_{6,\mathfrak{t}}^{\prime}=\frac{1}{4}\big{(}6G_{2,\mathfrak{t}}G_{6,\mathfrak{t}}-5G_{4,\mathfrak{t}}^{2}-4G_{8,\mathfrak{t}}\big{)},$
$\displaystyle
G_{8,\mathfrak{t}}^{\prime}=\frac{1}{2}(4G_{2,\mathfrak{t}}G_{8,\mathfrak{t}}-G_{4,\mathfrak{t}}G_{6,\mathfrak{t}}).$

Using a procedure outlined in [4, Section 4], we may use these equations to
arrive at a differential equation $C_{(4,4)}$ satisfied by
$y=G_{2,\mathfrak{t}}(\tau)$:

$\displaystyle C_{(4,4)}\colon\ y^{(4)}-5yy^{(3)}-9yy^{\prime
2}+(6y^{\prime}+6y^{2})y^{\prime\prime}=0.$ (4.2)

######  Remark 4.2.

The differential equation (4.2), as we will see shortly, has Bureau symbol
P$1$, meaning the differential equation admits solutions having a Laurent
expansion with leading divergence ${\sim}(\tau-\tau_{0})^{-1}$. This equation
is known to possess the Painlevé property, and appears in [10, equation
(2.3)].

These differential equations are natural analogues of the Chazy equation
associated to the quasimodular form of ${\rm PSL}(2,\mathbb{Z})$ [8, 24], and
also the higher-order Chazy eqautions associated to Hecke groups [4]. The
above systematic procedure for constructing higher-order nonlinear
differential equations motivates the following proposition.

######  Proposition 4.3.

Each triangle group is associated to a nonlinear differential equation of
order $m_{2}$, which we shall denote by $C_{\mathfrak{t}}$. Further, each term
in this differential equation is of order $2m_{2}+2$.

###### Proof.

The Ramanujan identities in terms of the combinations $G_{2k,\mathfrak{t}}$
are a system of $m_{2}$ first order differential equations in $m_{2}$
variables. The proposition follows by differential elimination. ∎

###  4.2 Painlevé analysis

We now briefly recall the main steps of the generalized Painlevé analysis,
following the work of [9]. Consider a differential equation represented by

$\displaystyle M[y]=0.$ (4.3)

We first propose a local solution $y_{c}$ (around a simple pole singularity)
to the nonlinear differential equation in Frobenius form. One then linearises
the equation around the solution

$\displaystyle\frac{\text{d}missing}{\text{d}\epsilon}M[y_{c}+\epsilon
w]\Big{|}_{\epsilon=0}=0,$ (4.4)

and look for solutions $w=w(\tau)$ also in the Frobenius form. In this case,
as in the original Chazy equation and in the higher order generalizations in
[4], we find that the linearised solution has poles of higher order than the
originally proposed solution to the differential equation. At first glance
these so-called “negative resonances” appear to make the leading order
solution unstable to linear perturbations. However, as shown in [9, 12], it is
possible to perform the stability analysis in such a situation by writing the
solution as a functional Frobenius series:

$\displaystyle y=\sum_{i=0}^{\infty}y_{i}\chi^{i-a}.$

Here, $\chi=\chi(\tau)$ satisfies the Riccati equation of the form (see [4]
for more details):

$\displaystyle\chi^{\prime}=1-\frac{S}{2}\chi^{2},$

and $y_{i}$ are the undetermined coefficients in the Frobenius ansatz. We find
that leading order ansatz satisfies the differential equation provided $a$ is
an integer solution to an indicial equation.555Integrality of solutions to an
indicial equation is one of the criteria laid out in [9], and ensures single-
valuedness of solutions. Once a particular branch of the indicial equation is
chosen, we can solve for the coefficients $y_{i}$. We then propose the
following ansatz for the next to leading order solution:

$\displaystyle w=\sum_{i=0}^{\infty}w_{i}\chi^{i-b}.$

We then substitute this into the linearized equation in (4.4) and solve for
$b$. The phenomenon of “negative resonances” is when $b$ is a positive integer
larger than $a~{}$. This appears to destabilize the leading order solution as
it appears to have a more singular behaviour near $\chi=0$. However on setting
the coefficients of the leading order solution to be the $y_{i}$ that we have
solved for, we find that the coefficients of the higher order negative
resonance vanishes and this implies that the negative resonance does not
destabilize the solution of the differential equation. Thus, following the
methods of [9, 12] we can claim that the Chazy equation for the triangle
groups of type $(M,M,\infty)$ satisfies the Painlevé property and this leads
us to the following proposition.

######  Proposition 4.4.

For a triangle group of type $\mathfrak{t}=(M,M,\infty)$, the associated Chazy
equation $C_{\mathfrak{t}}$ possesses the Painlevé property. The Eisenstein
series $y=E_{2,\mathfrak{t}}(\tau)$ that satisfies both the differential
equations has the following leading order expansion:

$\displaystyle
y=\frac{M}{M-2}\left(-\frac{2}{\chi}+\frac{S}{3}\chi+\frac{2S^{2}}{45}\chi^{3}+\frac{11S^{3}}{945}\chi^{5}+\cdots\right)\\!.$

The indicial equation for the resonances has the roots
$b\in\\{2,\dots,m+1\\}$.

This extends the results obtained for the Hecke groups in [4] to all isosceles
triangle groups with a cusp. Before concluding, we briefly discuss the status
of the Chazy equations for triangle groups of type $(m_{1},m_{2},\infty)$ with
$m_{1}\neq m_{2}$, and the Maier equation presented in Theorem 3.7. In the
former case, some of the resonance numbers $b$ are negative rational numbers.
We have observed that triangle groups of type $(m,m+k,\infty)$ have $k$ such
negative rational resonances. For the Maier equation, the analysis is more
intricate [20].666On revisiting our previous work, it appears that the our
analysis of the Maier equation [4, equation (4.5)] using the methods of
Conte–Fordy–Pickering [9] are insufficient to prove the Painlevé property. The
conclusions of [4, Proposition 5.1] pertaining to the higher-order Chazy
equations $C_{m}$ [4, Proposition 4.1], however, remain unchanged. A more
comprehensive analysis of the Maier equation is in preparation. In both these
cases, the criteria of [9, 12] that test for the Painlevé property are not
straightforwardly applicable. Work on extending these techniques is currently
underway.

##  Appendix A Generalized Halphen system

In this section we collect a few results about the Halphen system and the
Fourier expansions of the automorphic forms $E_{2k,\mathfrak{t}}$ associated
to triangle group of type $\mathfrak{t}=(m_{1},m_{2},\infty)$. In terms of the
parameters $(a,b,c)$ introduced in (2.2), the explicit solution of the Halphen
system in (2.1) is given by (see [11, Theorem 3]):

$\displaystyle
t_{1,\mathfrak{t}}(\tau)=\frac{1}{\alpha_{\mathfrak{t}}}(a-1)zQ(a,b;z)\,{}_{2}F_{1}(1-a,b;1;z){}_{2}F_{1}(2-a,b;2;z),$
$\displaystyle
t_{2,\mathfrak{t}}(\tau)-t_{1,\mathfrak{t}}(\tau)=\frac{1}{\alpha_{\mathfrak{t}}}Q(a,b;z)\,{}_{2}F_{1}(1-a,b;1;z)^{2},$
$\displaystyle
t_{3,\mathfrak{t}}(\tau)-t_{1,\mathfrak{t}}(\tau)=\frac{1}{\alpha_{\mathfrak{t}}}zQ(a,b;z)\,{}_{2}F_{1}(1-a,b;1;z)^{2}.$
(A.1)

The function ${}_{2}F_{1}$ is the Gauss hypergeometric function while the
function $Q(a,b;z)$ is given by

$\displaystyle Q(a,b;z)=\frac{\text{i}\pi(1-b)}{2\sin(\pi b)\sin(\pi
a)}(1-z)^{b-a}.$

The parameter $z$ is related to the standard hauptmodul $J_{\mathfrak{t}}$ of
the triangle group:

$\displaystyle z=\frac{1}{1-J_{\mathfrak{t}}}.$

The solutions in [11] differ from those in (A.1) by an overall constant
$\alpha_{\mathfrak{t}}$ that we are free to choose. We choose this coefficient
such that the expansion for $t_{2,\mathfrak{t}}$ begins with unit coefficient;
this uniquely fixes the coefficient $\alpha_{\mathfrak{t}}$ to be

$\displaystyle\alpha_{\mathfrak{t}}=-\frac{m_{1}+m_{2}+m_{1}m_{2}}{8m_{1}m_{2}}\csc\left[\frac{\pi}{2}\left(1-\frac{1}{m_{1}}+\frac{1}{m_{2}}\right)\right]\sec\left[\frac{\pi(m_{1}+m_{2})}{2m_{1}m_{2}}\right]\\!.$

With this choice of normalisation, the quasi-automorphic weight-$2$ Eisenstein
series $E_{2,\mathfrak{t}}$ has the following behavior as
$\tau\rightarrow\text{i}\infty$:

$\displaystyle E_{2,\mathfrak{t}}(\tau)=1+O(q).$

##  Appendix B Chazy equations

In this appendix we give examples of nonlinear differential equations
naturally associated to triangle groups and solved by their quasiautomorphic
weight-$2$ Eisenstein series. These differential equations were discussed in
Proposition 4.4 and are arrived at by differential elimination on the system
of equations governing the $G_{2k,\mathfrak{t}}$ defined in (4.1). For a
triangle group of type $\mathfrak{t}=(m_{1},m_{2},\infty)$, they are denoted
$C_{\mathfrak{t}}\equiv C_{(m_{1},m_{2})}$. The following equations all
possess the Painlevé property:

$\displaystyle C_{(3,3)}\colon\ y^{(3)}-2yy^{\prime\prime}+3y^{\prime 2}=0,$
$\displaystyle C_{(4,4)}\colon\ y^{(4)}-5yy^{(3)}-9yy^{\prime
2}+\big{(}6y^{\prime}+6y^{2}\big{)}y^{\prime\prime}=0,$ $\displaystyle
C_{(5,5)}\colon\ 75y^{(5)}-675yy^{(4)}+520y^{\prime\prime 2}-366y^{\prime
3}+2727y^{2}y^{\prime 2}$ $\displaystyle\phantom{C_{(5,5)}\colon\
}{}+y^{(3)}\big{(}220y^{\prime}+1959y^{2}\big{)}+\big{(}{-}2664yy^{\prime}-1818y^{3}\big{)}y^{\prime\prime}=0,$
$\displaystyle C_{(6,6)}\colon\ 3y^{(6)}-42y^{(5)}y-228yy^{\prime\prime
2}-552y^{3}y^{\prime 2}+228yy^{\prime
3}+y^{(4)}\big{(}214y^{2}-12y^{\prime}\big{)}$
$\displaystyle\phantom{C_{(6,6)}\colon\
}{}+\big{(}600y^{2}y^{\prime}-114y^{\prime
2}+368y^{4}\big{)}y^{\prime\prime}+y^{(3)}\big{(}57y^{\prime\prime}-34yy^{\prime}-468y^{3}\big{)}=0.$

### Acknowledgments

We thank Robert Conte for helpful correspondence, and Hossein Movasati for
valuable discussions. We would also like to thank the anonymous referees for
valuable comments and feedback. DPJ and MR are grateful to IMSc, Chennai for
hospitality. MR acknowledges support from the Infosys Endowment for Research
into the Quantum Structure of Spacetime. DPJ acknowledges support from SERB
grant CRG/2018/002835.

## References

  * [1] Ablowitz M.J., Chakravarty S., Hahn H., Integrable systems and modular forms of level 2, J. Phys. A: Math. Gen. 39 (2006), 15341–15353, arXiv:math.NT/0609210.
  * [2] Ablowitz M.J., Ramani A., Segur H., A connection between nonlinear evolution equations and ordinary differential equations of $P$-type. II, J. Math. Phys. 21 (1980), 1006–1015.
  * [3] Ashok S.K., Dell’Aquila E., Lerda A., Raman M., S-duality, triangle groups and modular anomalies in ${\mathcal{N}}=2$ SQCD, J. High Energy Phys. 2016 (2016), no. 4, 118, 52 pages, arXiv:1601.01827.
  * [4] Ashok S.K., Jatkar D.P., Raman M., Aspects of Hecke symmetry: anomalies, curves, and Chazy equations, SIGMA 16 (2020), 001, 26 pages, arXiv:1810.07919.
  * [5] Başar G., Dunne G.V., Ünsal M., Quantum geometry of resurgent perturbative/nonperturbative relations, J. High Energy Phys. 2017 (2017), no. 5, 087, 56 pages, arXiv:1701.06572.
  * [6] Bureau F.J., Sur des systèmes différentiels non linéaires du troisième ordre et les équations différentielles non linéaires associées, Acad. Roy. Belg. Bull. Cl. Sci. (5) 73 (1987), 335–353.
  * [7] Cheng M.C.N., Dabholkar A., Borcherds–Kac–Moody symmetry of ${\mathcal{N}}=4$ dyons, Commun. Number Theory Phys. 3 (2009), 59–110.
  * [8] Clarkson P.A., Olver P.J., Symmetry and the Chazy equation, J. Differential Equations 124 (1996), 225–246.
  * [9] Conte R., Fordy A.P., Pickering A., A perturbative Painlevé approach to nonlinear differential equations, Phys. D 69 (1993), 33–58.
  * [10] Cosgrove C.M., Higher-order Painlevé equations in the polynomial class. II. Bureau symbol $P1$, Stud. Appl. Math. 116 (2006), 321–413.
  * [11] Doran C.F., Gannon T., Movasati H., Shokri K.M., Automorphic forms for triangle groups, Commun. Number Theory Phys. 7 (2013), 689–737, arXiv:1307.4372.
  * [12] Fordy A., Pickering A., Analysing negative resonances in the Painlevé test, Phys. Lett. A 160 (1991), 347–354.
  * [13] Huber T., Differential equations for cubic theta functions, Int. J. Number Theory 7 (2011), 1945–1957.
  * [14] Jatkar D.P., Sen A., Dyon spectrum in CHL models, J. High Energy Phys. 2006 (2006), no. 4, 018, 32 pages, arXiv:hep-th/0510147.
  * [15] Maier R.S., Nonlinear differential equations satisfied by certain classical modular forms, Manuscripta Math. 134 (2011), 1–42, arXiv:0807.1081.
  * [16] Matsuda K., Differential equations involving cubic theta functions and Eisenstein series, Osaka J. Math. 57 (2020), 521–542, arXiv:1609.07481.
  * [17] Milnor J., On the $3$-dimensional Brieskorn manifolds $M(p,q,r)$, in Knots, Groups, and 3-Manifolds (Papers dedicated to the memory of R.H. Fox), Ann. of Math. Studies, Vol. 84, Princeton University Press, 1975, 175–225.
  * [18] Nehari Z., Conformal mapping, Courier Corporation, 2012.
  * [19] Persson D., Volpato R., Fricke S-duality in CHL models, J. High Energy Phys. 2015 (2015), 156, 55 pages, arXiv:1504.07260.
  * [20] Pikering A., Painlevé hierarchies and the Painlevé test, Theoret. and Math. Phys. 137 (2003), 1733–1742.
  * [21] Ramamani V., On some algebraic identities connected with Ramanujan’s work, in Ramanujan International Symposium on Analysis (Pune, 1987), Macmillan of India, New Delhi, 1989, 277–291.
  * [22] Raman M., Aspects of Hecke symmetry I: Ramanujan identities and inversion formulas, arXiv:1803.10224.
  * [23] Raman M., Bala Subramanian P.N., Chebyshev wells: periods, deformations, and resurgence, Phys. Rev. D 101 (2020), 126014, 18 pages, arXiv:2002.01794.
  * [24] Takhtajan L.A., A simple example of modular forms as tau-functions for integrable equations, Theoret. and Math. Phys. 93 (1992), 1308–1317.
  * [25] Tsanov V.V., Triangle groups, automorphic forms, and torus knots, Enseign. Math. 59 (2013), 73–113, arXiv:1011.0461.
  * [26] Zagier D., Quantum modular forms, in Quanta of Maths, Clay Math. Proc., Vol. 11, Amer. Math. Soc., Providence, RI, 2010, 659–675.