# A Study of Variations in Correlation Between Rotation Residual and
Meridional Velocity of Sunspot Groups

J. Javaraiah Bikasipura, Bengaluru-560 111, India.
Formerly working at Indian Institute of Astrophysics, Bengaluru-560 034,
India.
email: jajj55<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

We analyzed the combined 142 years sunspot-group data from Greenwich
Photoheliograpic Results (GPR) and Debrecen Photoheliographic Data (DPD) and
determined the yearly mean residual rotation rate and the meridional velocity
of sunspot groups in different $5^{\circ}$ latitude intervals. We find that
there exists a considerable latitude–time dependence in both the residual
rotation and the meridional motion. The residual rotation rate is found to be
$-120$ m $s^{-1}$ to $80$ m $s^{-1}$. In a large number of solar cycles the
rotation is some extent weaker during maxima than that of during minima. There
exist alternate bands of equatorward and poleward meridional motions. The
equatorward motion is dominant mostly around the maxima of solar cycles with
velocity 8–12 m s-1, whereas the poleward motion is dominant mostly around the
minima but with a relatively weak velocity, only 4–6 m s-1. The analysis of
the data during Solar Cycles 12–24 that are folded according to the years from
their respective epochs of maxima suggests the existence of equatorward
migrating alternate bands of slower and faster than average rotation within
the activity belt. This analysis suggests no clear equatorward or poleward
migrating bands of meridional motions. A statistically significant
anticorrelation exists between the meridional motion and residual rotation.
The corresponding linear-least-squares best-fit is found to be reasonably good
(slope, $-0.028\pm 0.008$, is about 3.5 times larger than its standard
deviation). As per the sign convention used for the meridional motions, the
significant negative value of the slope indicates the existence of a strong
angular momentum transport toward equator. The cross-correlation between the
slopes determined from the data in 3-year moving time intervals and yearly
mean sunspot number (SN) suggests that the slope leads SN by about 4 and 9
years. The Morlet wavelet spectrum of the slope suggests the existence of
$\approx$11-year periodicity in the slope almost throughout the data window,
but it was very weak during 1920–1940. We have also done cross-wavelet,
wavelet-coherence, and wavelet-phase difference analyrses of the slope and SN.
Overall the results suggest there exits a strong relationship between the
slope and amount of activity during a solar cycle. However, the correlation
between the cycle-to-cycle modulations in the slope and the amplitude of solar
cycle is found to be insignificant, indicating that there is no relationship
between the slope and strength of activity on a long-time scale (longer than
11-year period).

###### keywords:

Sun: Dynamo – Sun: surface magnetism – Sun: activity – Sun: sunspots

##  1 Introduction

Studies of variations in solar activity are important for a better
understanding the basic mechanism of solar activity and solar cycles of
various periods to make the long-term forecasts of space-weather and also may
be Earth’s climate (hath15). It is well believed that the solar dynamo
processes are responsible for solar activity and cycle. The dynamo processes
involve the formation of toroidal magnetic field due to the shearing of the
poloidal magnetic field by differential rotation and conversion of the
toroidal field into poloidal field of opposite polarity over the course of
approximately 11-year (dg06). The solar meridional circulation may play a
vital role to transfer the angular momentum and magnetic flux across the solar
latitudes and even can maintain the observed differential rotation (e.g.
schr85, and references therein). Therefore, the studies of the correlation
between the variations in the solar meridional flow and differential rotation
are important for understanding role of meridional flow in the variations of
solar differential rotationi, and hence the variations in solar activity. The
study of correlation between the meridional and rotational flows was started
by ward65, who analyzed Greenwich sunspot-group data during the period
1935–1944 and found a significant correlation between the angular and
meridional velocities of sunspot groups and interpreted it as meriodional
flows transfer angular momentum toward equator. Later many scientists studied
this correlation by using various data and methods. For example, pat91 from
the Greenwich sunspot data during the period 1874–1976 and Howard (1991; 1996)
using Mt. Wilson measurements of sunspot groups (1917–1985) and of plages
(1967–1985) confirmed the existence of a strong correlation between the
meridional and rotation velocities of these solar phenomena. However, some
authors found the existence of only small covariance of the meridional and
angular velocities derived from the data of different tracers (nesme93;
komm94; meun97). Theoretical models (e.g. gilm86) predicted the observed
correlation between the surface latitudinal and longitudinal motions as
reflection of equatorial angular momentum transport caused by Reynolds
stresses near the surface. There were criticisms of this theoretical
explanation on the basis that the observed motions may result from the well-
known expansion and contraction of sunspots along the tilted magnetic axes of
the sunspot groups. gilm84 argued that because the effect could be observed
for whole sunspot groups, at least some fraction of the observed correlation
must be due to Reynolds stresses near the solar surface and that the amount
was sufficient to account for the angular momentum transport required to
maintain the solar differential rotation. So far most of the authors
attributed the observed correlation between meridional and rotational motions
to the action of Reynolds stresses that demonstrated the presence of a net
transport of angular momentum towards the equator able to maintain the
differential rotation. Recently, sudar14 by analyzing Greenwich Photographic
Result (GPR) and Solar Observing Optical Network (SOON) sunspot-group data
covering the period from 1878 until 2011 and sudar17 by analyzing the Debrecen
Photoheliographic Data (DPD) during the period 1974–2016, found a
statistically significant correlation between sunspot groups’ meridional
velocities and rotation velocity residuals confirming the transfer of angular
momentum towards the equator.

Solar meridional flows vary during the solar cycle (snod87; komm93a; meun05;
svanda08; hath10). Meridional flows also seem to have a vital role in the
cause of solar 11-year period torsional oscillations (see snod92) and in solar
flux transport dynamo process (dg06). ju06 analyzed a large set of sunspot
group data (1874–2004) and found the existence of correlation (good in the
northern hemisphere and weak in the southern hemisphere) between the mean
solar-cycle variations of meridional flow and the latitude gradient term of
solar rotation. In the present study we analyze the combined updated sunspot-
group data reported in GPR during the period 1874–1976 and in DPD during the
period 1977–2017 and study variation in the correlation between the residual
rotation and meridional motion of sunspot groups and its relationship with
solar activity through cross-correlation and continuous- and cross-wavelet
analyses.

In the next section we describe the data analysis, in Section 3 we show the
results, and in Section 4 we summarize the conclusions and discuss them
briefly.

##  2 Data analysis

We have downloaded the daily sunspot-group data reported in GPR during the
period April 1874 – December 1976 and DPD during the period January 1977 –
June 2017 from the website fenyi.solarobs.unideb.hu/pub/DPD/. The details
about these data can be found in gyr10, bara16, and gyr17. These data contain,
beside other parameters, the date and time of observation, heliographic
latitude ($\lambda$) and longitude ($L$), and central meridian distance (CMD)
of a sunspot group for each day during its appearance on the solar disk. The
solar sidereal angular velocity $\omega$ (in degree day-1) and meridional
velocity $v_{mer}$ (in degree day-1) of a sunspot group are calculated by
using the latitudes and longitudes of the sunspot group measured at times
$t_{i}$ and $t_{i-1}$ during the life time (disk passage) of the sunspot group
as follows (here $t$ is the date $+$ fraction of the day corresponding to the
time of observation):
$\omega(\theta)=\frac{L_{i}-L_{i-1}}{t_{i}-t_{i-1}}+14^{\circ}.18$ and
$v_{mer}(\theta)=\frac{\lambda_{i}-\lambda_{i-1}}{t_{i}-t_{i-1}}$, where
$14^{\circ}.18$ day-1 is the Carrington rigid body rotation. In all our
earlier analyses we have assigned the velocity value to the mean of
$\lambda_{i-1}$ and $\lambda_{i}$. Following the suggestion by ok05 in jj20
and here we assigned the velocity value to the $\theta=\lambda_{i-1}$ (also
see sudar14). Each disk passage of a recurrent sunspot group is treated as an
independent sunspot group. Hence, we have considered all the sunspot groups
that have life times 2–12 days. We have not used the data on the days when the
$|{\rm CMD}|>75^{\circ}$. This reduces the foreshortening effect if any. In
addition, the data correspond to the absolute latitudinal drifts $>2^{\circ}$
day-1 and absolute longitudinal drifts $>3^{\circ}$ day-1 are excluded. This
reduces considerably the uncertainty in the derived results (ward65; jg95). In
the case of meridional velocity here we have used the following sign
convention: in both the northern and southern hemispheres positive and
negative values indicate the poleward and equatorward motions, respectively.
We have converted the meridional velocity that measured in degree day-1 into m
s-1 ($0.01\ {\rm degree\ day}^{-1}=1.4\ {\rm m\ s}^{-1}$, cf. rf91).

The data of sunspot groups in the northern and southern hemispheres are
combined. We binned the daily values of $\omega$ and $v_{mer}$ during each of
the years 1874–2017 into different $5^{\circ}$ latitude intervals within
$0^{\circ}$–$40^{\circ}$ sunspot latitude (absolute) belt. We determined the
mean values $\langle\omega(\theta)\rangle$ and $\langle
v_{mer}(\theta)\rangle$ of the daily values of $\omega$ and $v_{mer}$,
respectively, in each $5^{\circ}$ latitude interval during each year (note
that $\langle\ \rangle$ implies the mean over a time interval and $\theta$
represents the middle value of a latitude interval). In several latitudes
intervals, particularly in a large extent during the years corresponding to
the minima of the solar cycles, the numbers of velocity values are found to be
zero. In some latitude intervals they are found to be just one. The data in
all such latitude intervals of all the years are excluded. We determined mean
and standard error ($\sigma_{mer}$) of $\langle v_{mer}(\theta)\rangle$ of the
remaining all latitude intervals and in all years. The
$\langle\omega(\theta)\rangle$ and $\langle v_{mer}(\theta)\rangle$ in the
latitude intervals which correspond to $\langle
v_{mer}(\theta)\rangle>2.5\sigma$ are excluded. We fitted the values of
$\langle\omega(\theta)\rangle$ in the remaining all latitudes intervals during
all the years 1874–2017 into the standard law of differential rotation:
$\omega(\theta)=A+B\sin^{2}{\theta}$ degree day-1, where $A$ and $B$ represent
the equatorial rotation rate and latitude gradient of rotation, respectively,
and $\theta$ is the middle value of a 5∘ latitude interval. We obtained the
residual ($\Delta\langle\omega(\theta)\rangle$) of the mean angular velocity
$\langle\omega(\theta)\rangle$ in each latitude interval (i.e. at each value
of $\theta$) by subtracting the mean value of $\langle\omega(\theta)\rangle$
of that latitude interval from the value of $\langle\omega(\theta)\rangle$ of
the same latitude interval deduced from the law of the differential rotation.
We determined the correlation between $\langle v_{mer}(\theta)\rangle$ and
$\Delta\langle\omega(\theta)\rangle$ of all latitude intervals and all years.
The corresponding linear-least-squares best-fit is found to be mostly relevant
to the values of $\langle v_{mer}(\theta)\rangle\leq 20$ m s-1. This is
because the values of $\langle v_{mer}(\theta)\rangle$ that are $>20$ m s-1
are few (14 %). Hence, we have excluded the values of $\langle
v_{mer}(\theta)\rangle$ and $\langle\omega(\theta)\rangle$ in the latitude
intervals which are correspond to the $\langle v_{mer}(\theta)\rangle$$>20$ m
s-1. In addition, we have excluded the values of
$\langle\omega(\theta)\rangle$ and $\langle v_{mer}(\theta)\rangle$ that
correspond to $\Delta\langle\omega(\theta)\rangle>1^{\circ}$ day-1, which are
also very few (0.6 %). We repeated all the aforementioned calculations, i.e.
we fitted the remaining values of $\langle\omega(\theta)\rangle$ to the
standard law of differential rotation and obtained the differential rotation
law $\langle\omega(\theta)\rangle=(14.5\pm 0.01)-(2.2\pm
0.07)\sin^{2}{\theta}$ degree day-1. By using the values of
$\langle\omega(\theta)\rangle$ deduced from this law we obtained the values of
$\Delta\langle\omega(\theta)\rangle$ in all latitude intervals. We have
converted the values of the residual rotation
$\Delta\langle\omega(\theta)\rangle$ degree day-1 into $\Delta\langle
v_{rot}(\theta)\rangle$ m s-1 to have the same units for both $\Delta\langle
v_{rot}(\theta)\rangle$ and $\langle v_{mer}(\theta)\rangle$ for the sake of
convenience to the study correlation between these parameters. We determined
the correlation between $\langle v_{mer}(\theta)\rangle$ and $\Delta\langle
v_{rot}(\theta)\rangle$ and the corresponding linear-least-squares best-fit of
these parameters in all latitude intervals during the whole period 1874–2017
and also separately from the values in each of Solar Cycles 12–24. We have
used the epochs of minima (1878.958, 1890.204, 1902.042, 1913.623, 1923.623,
1933.707, 1944.124, 1954.288, 1964.791, 1976.206, 1986.707, 1996.624,
2008.958), as well as maxima (1883.958, 1894.042, 1906.123, 1917.623,
1928.290, 1937.288, 1947.371, 1958.204, 1968.874, 1979.958, 1989.874,
2001.874, 2014.288) and the maximum values of these cycles that were
determined by pesnell18 by using the 13-month smoothed monthly mean version-2
international sunspot number (SN) series. We determined the average behavior
of the correlation during Solar Cycles 12–24 by using the method of superposed
epoch analysis, i.e. in order to have a better statistics all the 13 solar
cycles data are superposed/combined according to the years from the respective
epochs of maxima of the solar cycles. Such analyses were done by many authors
(gilm84; balth86; jk99; jj03; ju06; braj06; sudar14). We have also binned the
values of $\Delta\langle v_{rot}(\theta)\rangle$ and $\langle
v_{mer}(\theta)\rangle$ into 3-year moving time intervals (3-year MTIs)
successively shifted by one year, namely 1874–1876, 1875–1877,…,2015–2017 and
determined the corresponding correlation and the linear-least-squares best-fit
in each of the 3-year MTIs.

Wavelet transform can be used to analyze time series that contain
nonstationary power at many different frequencies. That is, by wavelet
analysis it is possible to decompose a time series into time–frequency space.
Hence, one can determine both the dominant modes of variability and how those
modes vary in time. The cross-wavelet transform between two time series is
simply the product of the first complex wavelet transform with the complex
conjugate of the second. The cross-wavelet power spectrum can be used as a
quantified indication of the similarity of power between two time series. The
wavelet-coherency is a normalized time and scale (period) resolved measure for
the relationship between two time series. It is the square of the cross-
wavelet power spectrum normalized by the individual wavelet-power spectra.
This gives a quantity between 0 and 1. A value of 1 means the existence of
linear relationship between the two time series around a time on a scale. A
value of zero means vanishing correlation (no linear relationship).
Measurements of wavelet-phase difference between two time series yield
information on the phase delay between oscillations in the time series as a
function of frequency (for details see tc98). We have done Morlet wavelet
analysis for the time series of the slopes of the linear relationships in
3-year MTIs. We have also done the wavelet analysis for the yearly values of
the version-2 SN. The yearly average version-2 SN time series is downloaded
from www.sidc.be/silso/datafiles. The details of changes and corrections in
version-2 SN can be found in clette16. Similarities in the Morlet wavelet
spectra of the slope and SN are checked from cross-wavelet, wavelet-coherence,
and wavelet-phase difference analyses. We have used the IDL-codes of the
wavelet analyzes provided by tc98 and available at
paos.colorado.edu/research/wavelets.

Figure 1.: Variations in (a) yearly mean residual rotation rate
($\Delta\langle v_{rot}(\theta)\rangle$) and (b) yearly mean meridional
velocity ($\langle v_{mer}(\theta)\rangle$) in different $5^{\circ}$ latitude
intervals ($\theta$ is the middle value of a latitude interval) during
1874–2017. Northern and southern hemispheres’ data are folded. Positive and
negative values of meridional velocity represent poleward and equatorward
motions, respectively.

Figure 2.: Average variation in (a) yearly mean residual rotation rate
($\Delta\langle v_{rot}(\theta)\rangle$) and (b) in yearly mean meridional
velocity ($\langle v_{mer}(\theta)\rangle$) of sunspot groups in different
$5^{\circ}$ latitude intervals during different phase of Solar Cycles 12–24,
determined by superposing the yearly data (the values shown in Figure 1) of
the Solar Cycles 12–24 (Cycle 24 is incomplete) according to the years from
their respective epochs of maxima (the epochs of maxima are taken as zeros).
Northern and southern hemispheres’ data are folded. Positive and negative
values of meridional velocity represent poleward and equatorward motions,
respectively.

##  3 Results

###  3.1 Latitude-Time Dependence of $\Delta\langle v_{rot}\rangle$ and
$\langle v_{mer}\rangle$

Figure 1a shows the variations in yearly mean residual rotation
($\Delta\langle v_{rot}(\theta)\rangle$) whose values correspond to the
$\Delta\langle\omega(\theta)\rangle\leq 1^{\circ}$ day-1 ($\approx$ 140 m
s-1), and Figure 1b shows yearly mean meridional velocity ($\langle
v_{mer}(\theta)\rangle$) which has values only $\pm 20$ m s-1, of sunspot
groups in different $5^{\circ}$ latitude ($\theta$) intervals during
1874–2017.

As can be seen in Figure 1a, there is a considerable latitude–time dependence
in $\Delta\langle v_{rot}\rangle$. The values of $\Delta\langle
v_{rot}(\theta)\rangle$ are from $-120$ m $s^{-1}$ to $80$ m $s^{-1}$ (only a
few values are beyond this interval). In a large number of solar cycles the
rotation is some extent weaker during maxima than that of during minima. There
is an indication on high-to-low latitude migration in $\Delta\langle
v_{rot}(\theta)\rangle$ over 8–10-year period. Particularly, the band of
slower than average rotation seems to be migrating from $\approx 35^{\circ}$
latitude to around $5^{\circ}$ latitude during many solar cycles (not clearly
visible in the early solar cycles).

As can be seen in Figure 1b, there is also a considerable latitude–time
dependence in $\langle v_{mer}\rangle$. There seem to be alternate bands of
equatorial motion ($\langle v_{mer}(\theta)\rangle$ is negative) and poleward
motion ($\langle v_{mer}(\theta)\rangle$ is positive). Overall, the
equatorward motion is dominant with velocity 8–12 m s-1 mostly around maxima
of more solar cycles, whereas the poleward motion seems to be relatively weak
with velocity only 4–6 m s-1 and exists mostly around minima of more solar
cycles (also see ju06). However, in some solar cycles, e.g. 18 and 24, the
motion was seem to be equatorward during their whole periods. There is also
some indication of equatorward migration in the latitude–time dependence of
$\langle v_{mer}\rangle$ during a large number of solar cycles. The beginnings
and endings of the equatorward bands are not at exactly the beginnings and
endings of solar cycles. Poleward meridional flow with a speed of about 20 m
s-1 has been well established at the solar photospheric level (meun09).

Figures 2a and 2b show the average variations in yearly mean $\Delta\langle
v_{rot}\rangle$ and $\langle v_{mer}\rangle$ of sunspot groups in different
$5^{\circ}$ latitude intervals during Solar Cycles 12–24. This is determined
by superposing the yearly data (the values shown in Figure 1) of Solar Cycles
12–24 according to the years from their respective epochs of maxima. Northern
and southern hemispheres’ data are folded. Figure 2a shows the existence of
alternate bands of slower and faster than average rotation, i.e. alternate
bands of negative and positive values of $\Delta\langle
v_{rot}(\theta)\rangle$, within the activity belt. The about $10^{\circ}$ wide
slow band (residual rotational velocity is $-$70 to $-80$ ms-1) seems to be
originated around $35^{\circ}$ latitude and the narrow one (only about
$5^{\circ}$ wide) seems to be originated around $15^{\circ}$ latitude during
the minimum of a solar cycle and both migrated toward low latitudes. The
latter looks to be ended in the declining phase (before end) of solar cycle.
The speed of migration of former is much low up to middle of the declining
phases of solar cycles and then suddenly increased. It should be noted here
that there is a major drawback in a superposed epoch analysis. The lengths and
rise times of solar cycles are considerably different. The values at beyond
the epochs $-4$ and $4$ (cf. Figure 2a) represent for a few solar cycles
rather than an average of all solar cycles. In a similar analysis sudar14 did
not see any regularity in changes of pattern in $\Delta\langle
v_{rot}(\theta)\rangle$ over time.

As can be seen in Figure 2b there are bands of equatorward meridional motion
($\langle v_{mer}(\theta)\rangle$ is negative), separated by a band of
poleward motion in $26^{\circ}$–$29^{\circ}$ latitude interval. These bands
equatorward or poleward migrations are not clearly visible. The motions of
active regions largely represent the plasma motion in the sun’s subsurface
layers (rf96). During the whole decline phase and in the center of activity
belt the motion is mostly equatorward. Around the maxima of solar cycles at
the middle and low latitudes the motion is mostly poleward and in
$21^{\circ}$–$25^{\circ}$ and $30^{\circ}$–$35^{\circ}$ latitudes the motion
is equatorward. Close to the beginnings (epoch $-4$) of solar cycles, the
motion looks to be mostly poleward at all latitudes. Since the magnetic
structures of large/long-lived sunspot groups might anchor relatively deeper
than those of small sunspot groups (see jg97), the equatorward motions of
sunspot groups during maxima of solar cycles represent the motions in
relatively deeper layers of the sun than the poleward motions of sunspot
groups during minima of solar cycles and mostly in low latitudes. The overall
pattern of the motion is somewhat consistent with the concept of flux-
transport dynamo models (dg06). siva10 analyzed the Mt. Wilson and Kodaikanal
sunspot-group data and found the equatorward motion in all latitudes. These
authors have used only first days’ data of sunspot groups, whereas here we
have used the data in all days during the life times of the sunspot groups.
zhao04 by employing a time–distance technique of helioseismology found that
the residual meridional flows (the flows subtracted by the mean meridional
flow profile of 1996) converged toward the solar activity belts. Recently,
sudar14 noted that the meridional motion of sunspot groups is directed towards
the zone of solar activity. This result is largely confirmed here around the
maxima of solar cycles (see Figure 2b), however, the poleword motions (4–5
ms-1) are much weaker than the equatorword motions (8–10 ms-1). rf91 found
that sunspots groups tend to move away from the central latitude of activity.
This is not found here.

###  3.2 Comparison of Variations in $\Delta\langle v_{rot}(\theta)\rangle$
with Torsional Oscillations

The so-called torsional oscillation discovered by hl80 from Mt. Wilson
velocity data during the period 1967–1980 consists of alternating bands of
faster (or slower) than average rotation moving from high latitudes towards
the equator in $\approx$22-year time. The maximum amplitude of the torsional
oscillation is about 5 m s-1. The faster-than-average rotation band is located
on the equatorward side of the magnetic activity belt and a slower-than-
average rotation band is located on its poleward side. Here, obviously, both
the slower and faster than average rotation, i.e. bands of negative and
positive values of $\Delta\langle v_{rot}(\theta)\rangle$, exist within the
activity belt. The torsional oscillation pattern was also found from full-disk
magnetograms but with the maximum of amplitude of the magnetic pattern is
twice as large as that of the Doppler pattern (snod91; komm93b). The maximum
absolute value of $\Delta\langle v_{rot}(\theta)\rangle$, that is determined
here from sunspot-group data, is much larger than the amplitude of the
velocity as well as the magnetic torsional oscillations. It is well-known that
magnetic structures of sunspots are rotating 2–3 % faster than the surrounding
plasma (see jk02). This is commonly interpreted as the magnetic structures
anchored at different depths in the convective zone, being coupled to layers
rotating at a different velocity (see jg97). We know from helioseismology that
there is a rotation gradient just below the surface (antia02).

It has been claimed that the torsional pattern is present in sunspot rotation
(gm82; ttk83). ttk83 analyzed longitudinal and latitudinal motions of
recurrent sunspot groups using Greenwich data during 1874–1976. Although their
time resolution was coarse, they found some evidence for the 11-yr oscillation
in the sunspot zones with an amplitude of a few m s-1. gilm84 and balth86
observed faster and slower bands, but with no clear migratory character.
tern90 found an evidence of equatorward moving bands of torsional oscillation
through a very careful study of the sunspot drawings made during Cycle 21 at
Catania Astrophysical Observatory. The latitude-time dependent pattern of
$\Delta\langle v_{rot}(\theta)\rangle$ that is seen in Figures 1a and 2a is
largely similar to that aforementioned pattern. tern90 used only the data of
old sunspot groups, i.e. for each sunspot group only the data collected from
the 4th day of observation until the last observation available have been
taken into account. Here we have used all the available data of sunspot groups
whose life times were 2–12 days.

Meunier et al. (1997) analyzed the rotation of photospheric faculae obtained
at Meudon throughout Cycle 19 (1954–1964) and found bands of faster and slower
rotation rates with an amplitude of a few meters per second similar to the
torsional oscillations. maka97 studied the long-term variations of the
differential rotation of the solar large-scale magnetic field using synoptic
H$\alpha$ maps in the latitude zone from $+45^{\rm o}$ to $-45^{\rm o}$ in the
period of 1915–1990. In each solar cycle, they found a band of faster or
slower than average rotation moving from high to low latitudes. The slow band
roughly corresponds to the location of magnetic activity, which is in a large
extent similar as the patterns in $\Delta\langle v_{rot}(\theta)\rangle$ are
seen here.

Figure 3.: Scatter plot of meridional velocity ($\langle
v_{mer}(\theta)\rangle$ from $-20$ to $+20$ m s-1) versus residual rotation
rate ($\Delta\langle v_{rot}(\theta)\rangle$ $\approx-140$ to $+140$ m s-1)
determined by the combined data of sunspot groups during Solar Cycles 12–24
($r=-0.1$, number of data points $N=1388$, Student’s $\tau=3.73$, $P=0.001$).
Positive values of meridional velocity indicate poleward motions and negative
values indicate equatorward motions. The continuous line (red) represents the
linear best-fit (${\rm slope}=-0.028\pm 0.008$)

.

###  3.3 Correlation Between $\langle v_{mer}(\theta)\rangle$ and
$\Delta\langle v_{rot}(\theta)\rangle$

Figure 3 shows the correlation between $\langle v_{mer}(\theta)\rangle$ and
$\Delta\langle v_{rot}(\theta)\rangle$ determined from the whole data gone in
Figures 1a and 1b (the combined data of sunspot groups during Solar Cycles
12–24). As can be seen in this figure a statistically significant
anticorrelation between $\langle v_{mer}(\theta)\rangle$ and $\Delta\langle
v_{rot}(\theta)\rangle$ exists (also see sudar14). The correlation coefficient
($r$) is $-0.1$, which is significant at the 99.9% confidence level ($N=1388$,
Students’ t: $\tau=3.73$, $P=0.001$). The corresponding linear-least-squares
fit is reasonably good in the sense that the slope ($-0.028\pm 0.008$) is
about 3.5 times larger than its standard deviation. Based on the sign
convention used for meridional velocities of sunspot groups (Sect. 2 above),
the negative value of the slope can be considered as a measure of the strength
of angular momentum transport toward equator.

Table 1.: The values of intercept ($C$) and slope ($D$) of the linear relationship between $\langle v_{mer}(\theta)\rangle$ and $\Delta\langle v_{rot}(\theta)\rangle$ at each epoch (year) determined from the yearly mean values of Solar Cycles 12–24 superposed according to the years from their respective epochs of maxima, i.e. determined from the values shown Figure 2. The corresponding values of the correlation coefficient ($r$), Students’ t ($\tau$), probability ($Prob.$), and number of data points ($N$) are also given. The corresponding values of all these parameters determined from the combined values of $\langle v_{mer}(\theta)\rangle$ and $\Delta\langle v_{rot}(\theta)\rangle$ of all the latitude intervals during the entire period, $-6$ to $8$, are given in the last row. Epoch | $C$ | $D$ | $r$ | $\tau$ | $Prob.$ | $N$
---|---|---|---|---|---|---
$-6$ | $-0.67\pm 5.13$ | $-0.32\pm 0.18$ | $-0.87$ | 1.80 | 0.838 | 5
$-5$ | $0.74\pm 1.69$ | $0.06\pm 0.04$ | $0.26$ | 1.30 | 0.896 | 27
$-4$ | $1.93\pm 1.04$ | $-0.06\pm 0.03$ | $-0.26$ | 2.13 | 0.982 | 69
$-3$ | $-0.83\pm 0.77$ | $0.00\pm 0.02$ | $0.02$ | 0.18 | 0.573 | 124
$-2$ | $1.33\pm 0.70$ | $-0.01\pm 0.02$ | $-0.03$ | 0.37 | 0.643 | 148
$-1$ | $0.19\pm 0.56$ | $-0.04\pm 0.02$ | $-0.15$ | 1.90 | 0.970 | 169
0 | $-0.09\pm 0.52$ | $-0.07\pm 0.02$ | $-0.22$ | 2.98 | 0.998 | 172
1 | $0.72\pm 0.54$ | $0.00\pm 0.03$ | $-0.01$ | 0.07 | 0.530 | 161
2 | $-0.75\pm 0.58$ | $-0.13\pm 0.03$ | $-0.39$ | 5.05 | 1.000 | 149
3 | $-0.76\pm 0.63$ | $-0.05\pm 0.03$ | $-0.13$ | 1.46 | 0.927 | 127
4 | $-1.11\pm 0.73$ | $-0.08\pm 0.03$ | $-0.26$ | 2.86 | 0.997 | 113
5 | $0.32\pm 0.83$ | $-0.04\pm 0.04$ | $-0.12$ | 1.06 | 0.855 | 84
6 | $0.64\pm 1.01$ | $0.05\pm 0.03$ | $0.23$ | 1.95 | 0.972 | 70
7 | $1.83\pm 2.39$ | $0.15\pm 0.07$ | $0.44$ | 1.98 | 0.967 | 20
Whole | $-0.27\pm 0.74$ | $-0.085\pm 0.028$ | $-0.40$ | 3.06 | 0.998 | 50

We also determined correlation between $\langle v_{mer}(\theta)\rangle$ and
$\Delta\langle v_{rot}(\theta)\rangle$ at each epoch (year) in Figures 2a and
2b, i.e. during the average solar cycle of the superposed data of Solar Cycles
12–24. We have done linear-least-squares fit to the corresponding data. In
Table 1 we have given the values of the corresponding parameters, namely
correlation coefficient ($r$), the slope of best-fit linear relationship, etc.
In many years the value of $r$ found to be significant at above 95% confidence
level ($Prob.\geq 0.95$, i.e., $P\leq 0.05$), obtained from Student’s t-test.
The correlation of the combined data of all epochs is also found to be
statistically significant (see the last row in Table 1). Figure 4 shows the
variation in the slope during the average solar cycle. If we exclude the point
at the epoch $-6$, a 11-year period cycle pattern (anticorrelation with solar
cycle) seems to present in the slope, suggesting that there exists a
relationship between the slope and activity during solar cycles (also see
sudar14). However, the positive value of the slope that corresponds to the
minimum epochs of a few long cycles have a large uncertainty ($\sigma$,
standard deviation)

Figure 4.: Mean variation of the slope of the linear best-fit of meridional
velocity ($\langle v_{mer}(\theta)\rangle$) and residual rotation
($\Delta\langle v_{rot}(\theta)\rangle$) during Solar Cycles 12–24 determined
from the superposed epoch analysis (used the values of the slope given Table
1). Zero represents the maximum epoch of a solar cycle.

Figure 5a shows the variation in the slopes of the linear relationships of
$\langle v_{mer}(\theta)\rangle$ and $\Delta\langle v_{rot}(\theta)\rangle$
determined from the data in 3-year MTIs during the period 1874–2017. As in
some of our earlier analyses (e.g. jk99) we have revised the time series of
the slopes in 3-year MTIs by replacing those values of the slopes having
$\sigma$ greater than 1.7 times the median-$\sigma$ with the mean of their
respective neighbor values. They are at 1878, 1896, 1907, 1908, 1932, 1952,
1974, 2005, 2006, 2008, 2009, and 2010. (No value of $\sigma$ is found to be
greater than 2.5 times the median-$\sigma$.) In Figure 5a we have also shown
the variation in yearly mean values of SN during 1875–2016 (used the file
SN_y_tot_v2.0.txt that was downloaded from www.sidc.be/silso/datafiles). In
Figure 5b we have shown the cross-correlation between the slope and SN. As can
be seen in Figure 5a, there exist variations of the order of 11 years. In
fact, we get a reasonably good anticorrelation ($r=-0.27$, Student’s $t=3.28$,
$p=0.001$) between the slope (revised data) and SN, suggesting that there
exists a relationship between the slope and activity during solar cycles. The
3–4 cycles (30–40 years) periodic variations are relatively strong during the
times of early solar cycles. The existence of strong $\approx$11-year period
variation in the cross-correlation coefficient (see Figure 5b) strongly
indicates the existence of a strong relationship between the slope and solar
cycle. The peak at $lag=4$ implies that the slope leads SN by about four
years. There is also a relatively large peak at $lag=9$ (largest negative
value of the cross-correlation coefficient) suggesting that the slope leads SN
by about nine years. The peak at $lag=4$ suggests that the slope around the
preceding minimum of a solar cycle may be related to the strength of activity
around the maximum of the same solar cycle. A large negative/positive (less
negative) slope at the minimum of a solar cycle probably indicates that the
solar cycle will have a small/large maximum. This is some extent similar to
the relationship between the strength of polar fields at minimum of a solar
cycle and the solar cycle maximum (sch78). The peak at $lag=9$ suggests that a
large negative/positive slope at the epoch after 1–2-years from the maximum of
a solar cycle ($n$) indicates a small/large amplitude for the next solar cycle
($n+1$). This is some extent similar to the relationship between the sum of
the areas of the sunspot groups in the equatorial latitudes just after about
one-year from the maximum epoch of Solar Cycle $n$ and the maximum of Solar
Cycle $n+1$ (jj07). We checked whether we can use the aforementioned
relationships of the slope and SN for predicting the maximum of Solar Cycle
25. It is found impossible. That is, we get a very large uncertainty in the
predicted value. The aforementioned relations of the slope and SN indicate
that large poleward flows of magnetic flux (of decaying active regions) during
the declining phase of a solar cycle may enhance the strength of polar fields
at the following minimum. Large equatorward flows (including down flows at
active regions) of magnetic flux may enhance the strength of activity in the
equatorial latitudes just after one-year from the maximum of the solar cycle.

Figure 5.: (a) Plot of the values (cross) of the slopes of the linear best-
fits of the meridional velocity ($\langle v_{mer}(\theta)\rangle$) and
residual rotation ($\Delta\langle v_{rot}(\theta)\rangle$) determined from
sunspot group data in 3-year MTIs: 1874–1876, 1875–1877,$\dots$,2015–2017,
versus middle years of these intervals. The dashed curve (red) represents the
revised data after those values having $\sigma$ (standard deviation) greater
than 1.7$\times$the median $\sigma$ were replaced with mean of their
respective neighbor values. The continuous curve represents the variation in
yearly mean values of SN. The Waldmeier numbers of the solar cycles and the
corresponding epochs of the minima and the maxima of solar cycles are also
shown by the symbols $m$ and $M$, respectively. The horizontal continuous-line
represents the mean value over the whole period. (b) Plot of the coefficient
of cross-correlation (CC) between the slope and SN versus lag.

Figure 6.: (Upper panel) Wavelet power spectra (left) and global spectra
(right) of the slope of linear relationship between the meridional velocity
($\langle v_{mer}(\theta)\rangle$) and residual rotation ($\Delta\langle
v_{rot}(\theta)\rangle$) of sunspot groups during 1875–2016 shown in Figure
5a. (Lower panel) Wavelet power spectra (left) and global spectra (right) of
the annual mean SN during 1875–2016. The wavelet spectra are normalized by the
variances of the corresponding time series. The shadings are at the normalized
variances of 1.0, 3.0, 4.5, and 6.0. The dashed curves represent the 95 %
confidence levels deduced by assuming a white-noise process. The cross-hatched
regions indicate the cone of influence where edge effects become significant
(Torrence and Compo, 1998).

Figure 7.: Panels (a) and (b) show the cross-wavelet power spectrum and
wavelet-coherence spectrum, respectively, of sunspot number (SN) and the slope
of the linear relationship between meridional velocity ($\langle
v_{mer}(\theta)\rangle$) and residual rotation ($\Delta\langle
v_{rot}(\theta)\rangle$) during the period 1875 – 2016. The shadings are at
levels 1.0, 3.0, 4.5, and 6.0. The cross-hatched regions indicate the cone of
influence where edge effects become significant (Torrence and Compo, 1998).
Figure 8.: Wavelet phase difference between SN and the slope of the linear
relationship between meridional velocity ($\langle v_{mer}(\theta)\rangle$)
and residual rotation ($\Delta\langle v_{rot}(\theta)\rangle$) in 3-year MTIs
during the period 1875 – 2016. The cross-hatched regions indicate the cone of
influence where edge effects become significant (Torrence and Compo, 1998).

In Figure 6 we compare the Morlet wavelet spectra of the variations in the
slope and SN shown in Figure 5a. Obviously, 11-year periodicity is strongly
present in SN throughout the analyzed data window. This periodicity seems to
be present in the slope almost throughout the data window, but it was very
weak during 1920–1940. The 16–32-year periodicities were strongly present in
the slope during the time of early solar cycles. Figure 7 shows the cross-
wavelet and wavelet-coherence spectra of the slope and SN. The cross-wavelet
spectrum (Figure 7a) suggests that there exists a statistical significant
similarity in the temporal behaviors of the slope and SN only during the
period 1950–1970, because outside this interval the slope has slightly smaller
than 11-year periodicity (see Figure 6a). In Figure 7b there is a suggestion
on the existence of strong coherence between the $\approx$11-year variations
of the slope and SN before 1940 and 1980-onward. There is also a suggestion on
the existence of coherence between the 20–30-year periodicities of the slope
and SN from 1940-onward. There are episodes of 3–5-year periodicities in both
the slope and SN between 1880 and 2016. There is a suggestion of the coherence
in 40–50-year periodic variations in the slope and SN, but this signal is
within the area of cone-of-influence. That is, the signal of this periodicity
is not well resolved, hence this periodicity cannot be detected here due to
inadequate data. Figure 8 shows the spectra of the wavelet-phase difference of
the slope and SN. As can be seen in this figure, there exists
$\approx-180^{\circ}$ phase difference between the 11-year period variations
of the slope and SN before 1900 (the slope seems to be lead SN), during the
periods 1940–1960 and 1980–2000. A 20–30 year period variations of the slope
and SN seems to be having $\approx$ $-180^{\circ}$ phase difference throughout
1875–2016. This is consistent with the anticorrelation between the slope and
SN found above. There seem to be considerable phase differences also in the
aforementioned remaining all periodic variations. A $90^{\circ}$ difference
indicates that there is a considerable phase mixing in the corresponding
variations of the slope and SN.

###  3.4 Cycle-to-Cycle Variation in the Slope

In Table 2 we have given the values of the intercept and the slope of the
linear relationship between $\langle v_{mer}(\theta)\rangle$ and
$\Delta\langle v_{rot}(\theta)\rangle$ determined from the data during each of
Solar Cycles 12–24, and also from the combined data of all cycles. The
corresponding values of the correlation coefficient, Student’s t, and
probability are given. The cycle interval, number of data points ($N)$, and
maximum version-2 sunspot number ($R_{\rm M}$) are also given. In Figure 9 we
have shown the cycle-to-cycle variations in the slope during Solar Cycles
12–24. As can be seen in this figure, although in several solar cycles the
correlation is poor and the coefficients of linear best-fits have large
uncertainties (also see Table 2), there is a possibility of a considerable
variation in the slope on the time scale of about 3–4-cycles, suggesting the
existence of a 3–4-cycle periodicity in the slope. In Solar Cycles 15, 18, 19,
20, and 22 the slope has significant negative values. The value of Solar Cycle
24 is also some extent negative (this cycle is incomplete). In Solar Cycles
14, 17, and 21 the slope has positive values, but these values are not
significantly different from zero. The equatorward angular momentum transport
may be relatively much less (absent) in Solar Cycles 14, 17, and 21. The
correlation ($r=-0.28$) between the slope and the amplitude of solar cycle is
found to be insignificant, indicating that there is no relationship between
the slope and strength of activity on a long-time scale (longer than 11-year
period).

Table 2.: The values of intercept ($C$) and slope ($D$) of the linear relationship between $\langle v_{mer}(\theta)\rangle$ and $\Delta\langle v_{rot}(\theta)\rangle$ determined from the data during the whole period of each solar cycle and from the combined data of all Solar Cycles 12–24 (the last row). The corresponding values of the correlation coefficient ($r$), Students’ t ($\tau$) and probability ($Prob.$) are given. Cycle interval (Time), number of data points ($N$), maximum ($R_{\rm M}$) Version-2 sunspot numbers are also given. Cycle | Time | $R_{\rm M}$ | $C$ | $D$ | $r$ | $\tau$ | $Prob.$ | $N$
---|---|---|---|---|---|---|---|---
12 | 1878–1889 | 124.4 | $0.781\pm 0.819$ | $-0.023\pm 0.027$ | $-0.09$ | 0.85 | 0.802 | 101
13 | 1890–1901 | 146.5 | $-0.398\pm 0.859$ | $-0.028\pm 0.028$ | $-0.10$ | 0.97 | 0.834 | 107
14 | 1902–1912 | 107.1 | $0.075\pm 0.871$ | $0.029\pm 0.029$ | $0.11$ | 1.01 | 0.843 | 87
15 | 1913–1922 | 175.7 | $-1.422\pm 0.722$ | $-0.086\pm 0.030$ | $-0.29$ | 2.85 | 0.997 | 93
16 | 1923–1932 | 130.2 | $-0.381\pm 0.724$ | $-0.026\pm 0.024$ | $-0.11$ | 1.05 | 0.853 | 103
17 | 1933–1943 | 198.6 | $-0.388\pm 0.689$ | $0.006\pm 0.026$ | $0.02$ | 0.22 | 0.587 | 111
18 | 1944–1953 | 218.7 | $-1.000\pm 0.718$ | $-0.086\pm 0.029$ | $-0.27$ | 2.97 | 0.998 | 113
19 | 1954–1963 | 285.0 | $0.074\pm 0.698$ | $-0.047\pm 0.024$ | $-0.18$ | 1.99 | 0.976 | 119
20 | 1964–1975 | 156.6 | $1.064\pm 0.648$ | $-0.048\pm 0.023$ | $-0.18$ | 2.07 | 0.980 | 127
21 | 1976–1985 | 232.9 | $-0.436\pm 0.787$ | $0.018\pm 0.031$ | $0.05$ | 0.56 | 0.712 | 119
22 | 1986–1995 | 212.5 | $0.909\pm 0.724$ | $-0.050\pm 0.030$ | $-0.16$ | 1.69 | 0.953 | 113
23 | 1996–2007 | 180.3 | $1.383\pm 0.664$ | $-0.017\pm 0.030$ | $-0.05$ | 0.56 | 0.713 | 133
24 | 2008–2017 | 116.4 | $0.476\pm 0.949$ | $-0.053\pm 0.043$ | $-0.13$ | 1.23 | 0.889 | 88
All | 1878–2017 |  | $0.16\pm 0.20$ | $-0.028\pm 0.008$ | $-0.1$ | 3.73 | 0.9999 | 1388
Figure 9.: Cycle-to-cycle variation in the slope (open circle-dotted curve) of
the linear best-fit of meridional velocity ($\langle v_{mer}(\theta)\rangle$)
and $\Delta\langle v_{rot}(\theta)\rangle$ of sunspot groups during Solar
Cycles 12–24. The filled circle-continuous curve (red) represents the cycle-
to-cycle variation in solar-cycle amplitude ($R_{\rm M}$).

##  4 Conclusions and Discussion

We analyzed the combined 142 years GPR and DPD sunspot-group data and found
that there is a considerable latitude–time dependence in the yearly mean
residual rotation rate of sunspot groups. The yearly average residual rotation
rate is roughly from $-120$ m $s^{-1}$ to $80$ m $s^{-1}$. In a large number
of solar cycles the rotation is some extent weaker during maxima than that of
during minima. In many solar cycles there is an indication that bands of
residual rotation rate migrate from high to low latitudes over 8–10-year
periods. The analysis of the data during Solar Cycles 12–24 that are folded
according to the years from their respective epochs of maxima suggests the
existence of alternate bands of slower and faster than average rotation within
the activity belt. A $\approx 10^{\circ}$-wide slow band ($-$70 to $-80$ ms-1)
seems to be originated around $35^{\circ}$ latitude and a narrow one (only
$5^{\circ}$ wide) seems to be originated around $15^{\circ}$ latitude during
the minimum of a solar cycle and both migrated toward low latitudes. However,
their relationship with the well-known torsional oscillations is not clear.
The maximum absolute value of the residual rotation rate is much larger than
the amplitude of the velocity as well as the magnetic torsional oscillations.

There is also a considerable latitude–time dependence in the yearly mean
meridional motion of sunspot groups. There are alternate bands of equatorward
and poleward motions. Overall, the equatorward motion is dominant with
velocity 8–12 m s-1 mostly around maxima of more solar cycles, whereas the
poleward motion seems to be relatively weak with velocity only 4–6 m s-1
mostly around minima of more solar cycles. The analysis of the folded data of
Solar Cycles 12–24 suggests no clear equatorward or poleward migrating bands
of meridional motions.

There exists a statistically significant anticorrelation between the
meridional motion and residual rotation. The corresponding linear-least-
squares best-fit is reasonably good in the sense that the slope ($-0.028\pm
0.008$) is about 3.5 times larger than its standard deviation. As per the sign
convention used for the meridional velocities of sunspot groups, the
significant negative value of the slope can be considered as a measure of the
strength of angular momentum transport toward equator. This result is in a
qualitative agreement with previous results for sunspots (sudar14; sudar17).
However, our result has a lower amplitude, which might be a consequence of
different selection process. We can conclude that the results of present work
generally confirm earlier results (sudar14; sudar17) and represent, together
with them, a strong evidence that Reynolds stress is indeed the dominant
mechanism for maintaining solar differential rotation via angular momentum
transport towards equator.

The cross-correlation between the slopes determined from the data in 3-year
moving time intervals and yearly mean sunspot number (SN) suggests that the
slope leads SN by about 4 and 9 years. The Morlet wavelet spectrum of SN
suggests, obviously, 11-year periodicity is strongly present in SN throughout
the analyzed data window. The Morlet wavelet spectrum of the slope also
suggests the existence of $\approx$ 11-year periodicity in the slope almost
throughout the data window, but it was very weak during 1920–1940. The cross-
wavelet spectrum suggests that there exists a statistical significant
similarity in the temporal behaviors of the slope and SN only during the
period 1950–1970. The wavelet-coherence spectrum suggests that there exists a
strong coherence between the $\approx$ 11-year variations of the slope and SN
before 1940 and 1980-onward. There exists $\approx$ $-180^{\circ}$ phase
difference between the 11-year period variations of the slope and SN before
1900, during 1940–1960, and during 1980–2000. The 20–30-year period variations
of the slope and SN seem to be having $\approx$ $-180^{\circ}$ phase
difference throughout 1874–2017. However, this periodicity is very weak in
both the slope and SN. The above mentioned result that the slope leads SN by
about 4 and 9 years is consistent with the $180^{\circ}$ phase difference
between the $\approx$ 11-year variations and that between the $\approx$
20-year variations of the slope and SN. Overall these results suggest there
exits a strong relationship between the slope and amount of activity during a
solar cycle. However, the correlation between the cycle-to-cycle modulations
in the slope and the amplitude of solar cycle is found to be insignificant,
indicating that there is no relationship between the slope and strength of
activity on a long-time scale (longer than 11-year period).

Reynolds stress is produced by interaction of convective elements and the
Coriolis force which cause a correlation between longitudinal and latitudinal
velocity components. The motion of sunspots represents the motion of somewhat
deeper layer of the sun (see jj13). Recent local helioseismic studies suggest
that patterns of variations in the large scale flows are the same in both
quite and active regions and the flow patterns in active regions associated to
a deeper layers than those of quite regions (komm20). The anticorrelation
between the meridional motion and residual rotation of sunspot groups found
here implies the existence of couplings between these motions, somewhere in
the convection zone. Since Reynolds stress produced can maintain the
equatorial angular momentum transport and thus differential rotation (gilm86),
hence there exists correlation between the strength of activity and the slope
of linear relationship between the meridional motion and residual rotation
during a solar cycle. The existence of a $\approx 180^{\circ}$ phase
difference between the slope and activity during their 11-year and 20–30-year
periodic variations implies the existence of couplings between the strength of
magnetic activity and the angular momentum transport toward equator on these
time scales. However, our analysis suggest that there is no such relationship
on longer than these time scales.

##  5 Acknowledgments

The author thanks the anonymous reviewer for useful comments and suggestions.
The author acknowledges the work of all the people contribute and maintain the
GPR and DPD Sunspot databases. The sunspot-number data are provided by WDC-
SILSO, Royal Observatory of Belgium, Brussels. The wavelet software was
provided by C. Torrence and G. Compo and is available at
paos.colorado.edu/research/wavelets.

##  6 Conflict of interest

The authors declare that they have no conflicts of interest.

## References

  * Antia (2002) Antia, H.M.: 2002, in: Javaraiah, J. and Gokhale, M.H. (eds.) The Sun’s Rotation, Nova Science, New York, p. 81.
  * Balthasar, Vázquze, and Wöhl (1986) Balthasar, H., Vázquez, M., Wöhl, H.: 1986, Astron. Astrophys. 155, 87. ADS: http://adsabs.harvard.edu/abs/1986A%26A…155…87B
  * Baranyi et al. (2016) Baranyi, T., Győri, L., Ludmány, A.: 2016, Solar Phys. 291, 3081. DOI: 10.1007/s11207-016-0930-1
  * Brajša, Ruždjak, and Wöhl (2006) Brajša, R., Ruždjak, D., Wöhl, H.: 2006, Solar Phys. 237, 365. DOI: 10.1007/s11207-006-0076-7
  * Clette and Lefévre (2016) Clette, F., Lefévre, L.: 2016, Solar Phys. 291, 2629. DOI: 10.1007/s11207-016-1014-y
  * Dikpati and Gilman (2006) Dikpati, M., Gilman, P.A.: 2006, Astrophys. J. 649, 498. DOI: 10.1086/506314
  * Gilman (1986) Gilman, P.A.: 1986, In: Sturrock P.A., Holzer, T.E., Mihalas D.M., Ulrich R.K. (eds.), Physics of the Sun, 1, (Dordrecht: D. Reidal), P. 95.
  * Gilman and Howard (1984) Gilman, P.A., Howard, R.: 1984, Astrophys. J. 283, 385. DOI: 10.1086/162316
  * Godoli and Mazzuconi (1982) Godoli, G., Mazzuconi, F.: 1982, Astron. Astrophys. 116, 188.
  * Győri, Baranyi, and Ludmány (2010) Győri, L., Baranyi, T., Ludmány, A.: 2010, in Proc. Intern. Astron. Union 6, Sympo. S273, 2011, 403. DOI: 10.1017/s174392131101564X
  * Győri et al. (2017) Győri, L., Ludmány, A., Baranyi, T.: 2017, Mon. Not. Roy. Astron. Soc. 465, 1259. DOI: 10.1093/mnras/stw2667
  * Hathaway (2015) Hathaway, D.H.: 2015, Liv. Rev. Solar. Phys. 12(4), 1. DOI: 10.1007/1rsp-2015-4
  * Hathaway and Rightmire (2010) Hathaway, D.H., Rightmire, L.: 2010, Sci. 327, 1350. DOI: 10.1126/science.1181990
  * Howard (1991) Howard, R.F.: 1991, Solar Phys., 135, 327. DOI: 10.1007/BF00147504
  * Howard (1996) Howard, R.F.: 1996, Annu. Rev. Astron. Astrophys. 34, 75. DOI: 10.1146/annurev.astro.34.1.75
  * Howard and LaBonte (1980) Howard, R. LaBonte, B.J.: 1980, Astrophys. J. 239, L33. DOI: 10.1086/183286
  * Javaraiah (2007) Javaraiah, J.: 2007, Mon. Not. R. Astron. Soc. 377, L34. DOI: 10.1111/j.1745-3933.2007.00298.x
  * Javaraiah (2003) Javaraiah, J.: 2003, Solar Phys. 212, 23. DOI: 10.1023/A:1022912430585
  * Javaraiah (2013) Javaraiah, J.: 2013, Solar Phys., 287, 197.
  * Javaraiah (2020) Javaraiah, J.: 2020, Solar Phys., 295, 170. DOI: 10.1007/s11207-020-01740-x
  * Javaraiah and Gokhale (1995) Javaraiah, J., Gokhale, M.H.: 1995, Solar Phys. 158, 173. DOI: 10.1007/BF00680841
  * Javaraiah and Gokhale (1997) Javaraiah, J., Gokhale, M.H.: 1997, Astron. Astrophys. 327, 795.
  * Javaraiah and Komm (1999) Javaraiah, J., Komm, R.W.: 1999, Solar Phys. 184, 41. DOI: 10.1023/A:1005028128077
  * Javaraiah and Komm (2002) Javaraiah, J., Komm, R.W.: 2002, in: Javaraiah, J. and Gokhale, M.H. (eds.) The Sun’s Rotation, Nova Science, New York, p. 23.
  * Javaraiah and Ulrich (2006) Javaraiah, J., Ulrich, R.K.: 2006, Solar Phys. 237, 245. DOI: 10.1007/s11207-006-0130-5
  * Komm, Howard, and Harvey (1993a) Komm, R.W., Howard, R.F., Harvey, J.W.: 1993a, Solar Phys., 147, 207\. DOI: 10.1007/BF00690713
  * Komm, Howard, and Harvey (1993b) Komm, R.W., Howard, R.F., Harvey, J.W.: 1993b, Solar Phys. 143, 19. DOI: 10.1007/BF00619094
  * Komm, Howard, and Harvey (1994) Komm, R.W., Howard, R.F., Harvey, J.W.: 1994, Solar Phys. 151, 15. DOI: 10.1007/BF00654078
  * Komm, Howe, and Hill (2020) Komm, R., Howe, R., Hill, F.: 2020, Solar Phys. 295, 47. DOI: 10.1007/s11207-020-01611-5
  * Makarov, Tlatov, and Callebaut (1997) Makarov, V.I., Tlatov, A.G., Callebaut, D.K..:1997, Solar Phys. 170, 373. DOI: 10.1023/A.1004995826593
  * Meunier (2005) Meunier, N.: 2005, Astron. Astrophys. 442, 693. DOI: 10.1051/0004-6361:20053
  * Meunier and Zhao (2009) Meunier, N., Zhao, J.: 2009, Space Sci. Rev. 144, 127. DOI: 10.1007/s11214-008-9472-7
  * Meunier, Nesme-Ribes, and Collin (1997) Meunier, N., Nesme-Ribes, E., Collin, B.: 1997, Astron. Astrophys. 319, 683.
  * Nesme-Ribes, Ferreira, and Vince (1993) Nesme-Ribes, E., Ferreira, E.N., Vince, I.: 1993, Astron. Astrophys. 276, 211.
  * Olemskoy and Kitchatinov (2005) Olemskoy, S.V., Kitchatinov, L.L.: 2005, Astron. Lett. 31, 706. DOI: 10.1134/1.2075313
  * Paternó et al. (1991) Paternó, I., Spadaro, D., Zappalá, R.A., Zuccarello, F.: 1991, Astron. Astrophys., 252, 337.
  * Pesnell (2018) Pesnell, W.D.: 2018, Space Weather 16, 1997. DOI: 10.1029/2018SW002080
  * Schtten et al. (1978) Schatten, K.H., Scherrer, P.H., Svalgaard, L., Wilcox, J.M. 1978, Geophys. Res. Lett., 5, 411. DOI: 10.1029/GL005i005p00411
  * Schröter (1985) Schröter, E.H.: 1985, Solar Phys. 100, 141. DOI: 10.1007/BF00158426
  * Sivaraman et al. (2010) Sivaraman, K.R., Sivaraman, H., Gupta, S.S., Howard, R.F.: 2010, Solar Phys. 266, 247. DOI: 10.1007/s11207-010-9620-6
  * Snodgrass (1987) Snodgrass, H.B.: 1987, Astrophys. J. 316, L91. DOI: 10.1086/184899
  * Snodgrass (1991) Snodgrass, H.B.: 1991, Astrophys. J. 383, L85. DOI: 10.1086/186248
  * Snodgrass (1992) Snodgrass, H.B.: 1992, In: Harvey, K.L.(ed.) The Solar Cycle, CS-27, Astron. Soc. Pac., San Francisco, 205.
  * Sudar et al. (2014) Sudar, D., Skokić, I., Ruždjak, D., Brajsa, R., Wöhl, H.: 2014, Mon. Not. Roy. Astron. Soc. 439, 2377. DOI: 10.1093/mnras/stu099
  * Sudar et al. (2017) Sudar, D., Brajša, R., Skokić, I., Poljančć Beljan, I., Wöhl, H.: 2017, Solar Phys. 292, 86. DOI: 10.1007/s11207-017-1105-4
  * Švanda et al. (2008) Švanda, M., Klvaňa, M. Sobotka, V. Bumba: 2008, Astron. Astrophys. 477, 285. DOI: 10.1051/0004-6361:20077
  * Ternullo (1990) Ternullo, M.: 1990, Solar Phys. 127, 29. DOI: 10.1007/BF00158512
  * Torrance and Compo (1998) Torrence, Ch., Compo, G.P.: 1998, Bull. Am. Meteor. Soc. 79, 61. DOI: 10.1175/1520-0477(1998)079¡0061,APGTWA¿2.0.CO;2
  * Tuominen, Tuominen, and Kyöläinen (1983) Tuominen, J., Tuominen, I., Kyöläinen, J.: 1983, Mon. Not. Roy. Astron. Soc. 205, 691. DOI: 10.1093/mnras/205.3.691
  * Ward (1965) Ward, F.: 1965, Astrophys. J. 141, 534. DOI: 10.1086/148143
  * Zhao and Kosovichev (2004) Zhao, J., Kosovichev, A.G.: 2004, Astrophys. J. 603, 776. DOI: 10.1086/381489


# A Rapidly Converging Machin-like Formula for $\pi$

Oleg S. Alferov

(November 2023)

###### Abstract

We present a simple recurrent formula to generate the Machin-like expression
for calculating $\pi/4$. The method works for any denominator in the starting
term and always provides a finite decomposition. We show that the terms in the
Machin-like formula decrease so rapidly that the Lehmer’s measure can be made
arbitrarily small only by selecting the first term.

We introduce the concept of the partial Machin-like formula. While the growth
of the integer numbers may quickly render the computer implementation
impractical, the same reason restricts the total contribution of the high
terms. If the required precision is known in advance, the subset of the
expression may be selected to satisfy it.

We also present the Python program to compute the terms of the Machin-like
formula (full and partial), and its Lehmer’s measure.

Keywords: $\pi$, arctangent, finite series, Machin-like formula, Lehmer’s
measure.

##  1 Background

In the 1706, John Machin discovered the fact that the $\pi$ number can be
expressed as the finite sum of arctangents taken from certain fractions:

$\frac{\pi}{4}=4\arctan\frac{1}{5}-\arctan\frac{1}{239}.$

Since the Maclaurin series

$\arctan x=x-\frac{x^{3}}{3}+\frac{x^{5}}{5}-\frac{x^{7}}{7}+\dots$

converges fast if $x<1$ is small, Machin’s discovery made it possible to
calculate $\pi$ with high precision even before the advent of computers.

To honor his discovery, the identities in the form

$\frac{\pi}{4}=\sum_{k=0}^{N}m_{k}\arctan\frac{1}{q_{k}},\quad\hbox{where
$m_{k}\in\mathbb{Z}$ and $q_{k}\in\mathbb{N}$}$ (1)

are called the Machin-like formulas. Since then, many new identities have been
found. To estimate, in modern terms, the computational complexity of
expressing $\pi$ by the formula (1), Lehmer introduced the following measure
in 1938:

$\lambda=\sum_{k=0}^{N}\frac{1}{\log_{10}q_{k}}.$ (2)

The smaller $\lambda$, the easier is the computation of $\pi$. [1, 2, 3]

##  2 General facts and notations

Let us remind readers of a few well-known facts and simple statements. For
those who are curious, either the formal proof is provided in the Supplement,
or the reference is given to the source where the fact can be found.

###  2.1 Evolution of the arctangent argument

Let $\tan\alpha=a/b$ and $\tan\beta=c/d$ be rational numbers. Then
$\tan(\alpha\pm\beta)$ is also the rational number, immediately expressed via
trigonometric identities. Because of that, an arbitrary sum of arctangents
from rational numbers is an arctangent from the rational fraction again. [4]

$\begin{array}[]{l}\tan(\alpha\pm\beta)=\displaystyle\frac{\tan\alpha\pm\tan\beta}{1\mp\tan\alpha\tan\beta},\quad\Rightarrow\\\\[10.00002pt]
\arctan\displaystyle\frac{a}{b}\pm\arctan\displaystyle\frac{c}{d}=\arctan\displaystyle\frac{ad\pm
bc}{bd\mp ac}.\end{array}$ (3)

The shortest form of the Machin-like formula (1) provides the connection
between $\pi$ and arctangents:

$\frac{\pi}{4}=\arctan\frac{1}{1}.$ (4)

However, (4) is the only single-term Machin-like formula. Fact 1:

$\frac{\pi}{4}\neq m\arctan\displaystyle\frac{1}{q_{0}},\quad\hbox{$\forall$
$m$, $q_{0}\in\mathbb{N}$, if $q_{0}>1$.}$

We are going to subtract arctangents of the selected fractions from the right-
hand side of (4), one by one. Since $1\in\mathbb{Q}$, the difference will be
another arctangent of a rational number on each step:

$\displaystyle\frac{\pi}{4}=\arctan\displaystyle\frac{1}{1}=\left(\pm\arctan\displaystyle\frac{1}{q_{0}}+...\pm\arctan\displaystyle\frac{1}{q_{n-1}}\right)+\arctan\displaystyle\frac{A_{n}}{B_{n}}.$
(5)

The sum in the brackets represents the Machin-like formula that we are
building. (If some of $q$-s repeat, we can group them into one term.) The
rightmost term outside brackets is the “remainder”. If we add another term to
the sum, the “remainder” evolves:

$\begin{array}[]{c}\arctan\displaystyle\frac{A_{n}}{B_{n}}\mp\arctan\displaystyle\frac{1}{q_{n}}=\arctan\displaystyle\frac{A_{n+1}}{B_{n+1}},\\\\[11.99998pt]
\left\\{\begin{array}[]{l}A_{n+1}=q_{n}A_{n}\mp B_{n},\\\
B_{n+1}=q_{n}B_{n}\pm A_{n},\end{array}\right.\end{array}$ (6)

and our expression for extending the Machin-like formula becomes

$\displaystyle\frac{\pi}{4}=\left(\pm\arctan\displaystyle\frac{1}{q_{0}}+...\pm\arctan\displaystyle\frac{1}{q_{n-1}}\pm\arctan\displaystyle\frac{1}{q_{n}}\right)+\arctan\displaystyle\frac{A_{n+1}}{B_{n+1}}.$

We would also have to reduce the fraction by $\gcd(A_{n+1},B_{n+1})$, but it
is irrelevant to our study.

The series starts with values $A_{0}=1$ and $B_{0}=1$. On each step, we use
$q_{n}$, extend the Machin-like expression with the next term, and get the
“remainder” corresponding to the pair $A_{n+1}$, $B_{n+1}$. Our goal is to
select the sequence $\\{q_{n}\\}$ in a way that the “remainder” becomes $0$ in
the end, in other words, $A_{N+1}=0$. After that, the iterations stop, and we
get the new Machin-like formula.

###  2.2 Nearest integer

The floor function, the nearest integer less than or equal to the argument,
$y=\left\lfloor x\right\rfloor$ is well defined. In $\mathbb{Q}$, it
ultimately comes from Euclid’s division lemma [5]. (Let $x=a/b$, then
$\exists$ $\alpha,\beta\in\mathbb{Z}$, $0\leq\beta<b$, such as $a=\alpha
b+\beta$, then $\left\lfloor x\right\rfloor=\alpha$ and $0\leq x-\left\lfloor
x\right\rfloor=\beta/b<1$, Q.e.d.) In $\mathbb{R}$, the floor can be
introduced via approximations of $x$ in $\mathbb{Q}$.

The ceiling function, the nearest integer greater than or equal to the
argument, can be defined via the floor function: $\left\lceil
x\right\rceil=-\left\lfloor-x\right\rfloor$. Properties of the floor and
ceiling functions include: [6]

$\left\\{\begin{array}[]{ll}\left\lfloor x\right\rfloor=x=\left\lceil
x\right\rceil,&x\in\mathbb{Z},\\\\[1.99997pt] \left\lfloor
x\right\rfloor<x<\left\lceil x\right\rceil=\left\lfloor
x\right\rfloor+1,&x\notin\mathbb{Z}.\end{array}\right.$ (7)

Interestingly, the notation for the nearest integer is not so common. Let us
use the symbol $y=\left\lfloor x\right\rceil$ and the following intuitive
definition.

Definition. The nearest integer to $x\in\mathbb{R}$ is either $\left\lfloor
x\right\rfloor$, or $\left\lceil x\right\rceil$, whichever is closer to $x$
(also, exact halves are rounded up):

$\left\lfloor x\right\rceil=\left\\{\begin{array}[]{ll}\left\lfloor
x\right\rfloor,&\left|\left\lfloor x\right\rfloor-x\right|<\left|\left\lceil
x\right\rceil-x\right|,\\\\[1.99997pt] \left\lceil
x\right\rceil,&\hbox{otherwise.}\end{array}\right.$ (8)

Fact 2. $\left|\left\lfloor x\right\rceil-x\right|\leq 1/2$.

Corollary 1. If $x\in\mathbb{R}$ is given, and $y\in\mathbb{Z}$ is variable,
the minimum of $\left|y-x\right|$ equals to $\left|\left\lfloor
x\right\rceil-x\right|$. (For proof, substitute $y=\left\lfloor
x\right\rceil+v$.)

Corollary 2. Let $a,b\in R$ be given, $b>0$, and $m\in\mathbb{Z}$ is variable.
Then $m=\left\lfloor a/b\right\rceil$ provides the minimum for expression
$\left|a-mb\right|$.

###  2.3 Miscellaneous

For studying the partial Machin-like formulas, we need simplified means to
estimate the higher terms in expression (1).

Fact 3. Maclaurin series error for $x\in[0,1)$. See, for example, [7].

$\arctan
x=\sum_{k=0}^{K}(-1)^{k}\displaystyle\frac{k^{2k+1}}{2k+1}+\varepsilon,\quad\hbox{and}\quad\left|\varepsilon\right|\leq\displaystyle\frac{x^{2K+3}}{2K+3}.$
(9)

Corollary. $x-x^{3}/3<\arctan x<x$ for all $x\in(0,1)$.

Fact 4. Notable limit. See, for example, [8].

$\lim_{x\to\infty}\left(1+\displaystyle\frac{a}{x}\right)^{x}=e^{a}$

##  3 The idea

Let us informally present the method to build a parametric Machin-like formula
and postpone the strict proof until the next section. We select the starting
$q_{0}$ and use identity (6) to subtract the first term of (1) until the
difference is positive:

$\displaystyle\frac{\pi}{4}-m\arctan\displaystyle\frac{1}{q_{0}}=\arctan\displaystyle\frac{A_{1}}{B_{1}},\quad\hbox{such
as}\quad 0<\displaystyle\frac{A_{1}}{B_{1}}<\displaystyle\frac{1}{q_{0}}.$

After the first term is established, we extend the expression by adding the
next terms. We use the recurrent formula (6) again. Let us keep $A_{n}\geq 0$
$\forall\,n$, and require that the sequence $\\{A_{n}\\}$ is strictly
decreasing:

$0\leq A_{n+1}=q_{n}A_{n}-B_{n}<A_{n}.$

Given that $A_{n}>0$, it is equivalent to

$\displaystyle\frac{B_{n}}{A_{n}}\leq
q_{n}<\displaystyle\frac{B_{n}}{A_{n}}+1\quad\Rightarrow\quad
q_{n}=\left\lceil\displaystyle\frac{B_{n}}{A_{n}}\right\rceil.$

Notice that $q_{n}>0$ and $B_{n+1}>B_{n}>0$, therefore $q_{n+1}>q_{n}$. For
this reason, expression (1) has $m_{0}=m$, and $m_{k}=1$ for all $k>0$.

That is all. We calculate the next term based on the “remainder” until
$A_{n+1}=0$. Since the sequence $\\{A_{n}\\}$ of integers is limited from
below and is strictly decreasing, it is finite. Any starting $q_{0}>1$
generates a Machin-like identity in which terms can be made arbitrarily small.

Unfortunately, this method may have slow convergence, and Lehmer’s measure of
the resultant identity may be high. We need stronger constraints to improve
our identity.

##  4 The series

The series above uses all positive coefficients. Let us modify the approach
and allow negative terms in the series. We are now looking for the identity in
the form

$\displaystyle\frac{\pi}{4}=m\arctan\displaystyle\frac{1}{q_{0}}+\sum_{n=1}^{N}\delta_{n}\arctan\displaystyle\frac{1}{q_{n}},\quad\hbox{where
$q_{n}>0$, and $\delta_{n}=\pm 1$.}$ (10)

###  4.1 Starting term

Let $q_{0}>1$ be the first denominator. Let us rewrite (6) for the case of
subtracting the first arctangent term from $\pi/4$ multiple times:

$\displaystyle\frac{\pi}{4}-m\arctan\displaystyle\frac{1}{q_{0}}=\arctan\displaystyle\frac{a_{m}}{b_{m}}$
(11) $\left\\{\begin{array}[]{l}a_{0}=1\\\ b_{0}=1\\\
a_{k+1}=q_{0}a_{k}-b_{k}\\\ b_{k+1}=q_{0}b_{k}+a_{k}\\\ \end{array}\right.$
(12)

We need the absolute value of the “remainder” to be as small as possible.
According to the Fact 2.1, it will be nonzero. With fixed $q_{0}$, because of
Fact 2.2, Corollary 2, it happens at the nearest integer of the fraction

$\gamma=\left.\displaystyle\frac{\pi}{4}\right/\arctan\displaystyle\frac{1}{q_{0}},\quad
m=\left\lfloor\gamma\right\rceil\quad\Rightarrow\quad
0<\left|\arctan\displaystyle\frac{a_{m}}{b_{m}}\right|\leq\displaystyle\frac{1}{2}\arctan\displaystyle\frac{1}{q_{0}}\quad\hbox{
}$ (13)

For practical calculations, we cannot use real numbers or trigonometric
functions. Let us construct the method for computing $m$ using the recurrent
equations (12). Notice that the right-hand side of (11) eventually becomes
negative. Specifically, the “remainder” is still positive at
$m_{(-)}=\left\lfloor\gamma\right\rfloor$, and negative at
$m_{(+)}=\left\lceil\gamma\right\rceil$. Proof. By the floor and ceiling
properties (7),

$0<\gamma-m_{(-)}<1\quad\Rightarrow\quad
0<\displaystyle\frac{\pi}{4}-m_{(-)}\arctan\displaystyle\frac{1}{q_{0}}<\arctan\displaystyle\frac{1}{q_{0}}$
$m_{(+)}=m_{(-)}+1\quad\Rightarrow\quad\displaystyle\frac{\pi}{4}-m_{(+)}\arctan\displaystyle\frac{1}{q_{0}}<0,\quad\hbox{{\it
Q.e.d.}}$

Because the arctangent is monotonically increasing function,

$\arctan\displaystyle\frac{a_{m_{(-)}}}{b_{m_{(-)}}}>0,\quad\arctan\displaystyle\frac{a_{m_{(+)}}}{b_{m_{(+)}}}<0\quad\Rightarrow\quad\displaystyle\frac{a_{m_{(-)}}}{b_{m_{(-)}}}>0,\quad\displaystyle\frac{a_{m_{(+)}}}{b_{m_{(+)}}}<0$

Notice from equations (12), that when $a_{k+1}$ becomes negative for the first
time, all of $a_{k}$, $b_{k}$, and $b_{k+1}$ are still positive. In other
words, $m_{(-)}$ and $m_{(+)}=m_{(-)}+1$ are defined by the place where the
sequence $\\{a_{k}\\}$ changes sign:

$a_{m_{(-)}}>0,\quad a_{m_{(+)}}<0.$ (14)

Now, select between $m=m_{(-)}$ and $m=m_{(+)}$. The $(+)$ shall be selected
if the $(+)$ “remainder” is less than the $(-)$ one:

$\arctan\left|\displaystyle\frac{a_{m_{(-)}}}{b_{m_{(-)}}}\right|>\arctan\left|\displaystyle\frac{a_{m_{(+)}}}{b_{m_{(+)}}}\right|\quad\Rightarrow\quad
m=m_{(+)},$

or, taking into account (14), the same as

$m=\left\\{\begin{array}[]{ll}m_{(+)},&a_{m_{(-)}}b_{m_{(+)}}+a_{m_{(+)}}b_{m_{(-)}}>0,\\\\[1.99997pt]
m_{(-)},&\hbox{otherwise.}\end{array}\right.$ (15)

###  4.2 Reducing series

We are building the Machin-like expression in form (10), adding terms one by
one. To simplify tracking the terms and “remainders” in (10) with (6), we
introduce separate variable $\delta_{n}$ for the sign of the “remainder”. (The
arctangent function is odd.) Let us rewrite expression (5) for the “remainder”
on each step:

$\displaystyle\frac{\pi}{4}=m\arctan\displaystyle\frac{1}{q_{0}}+\sum_{k=1}^{n-1}\delta_{k}\arctan\displaystyle\frac{1}{q_{k}}+\delta_{n}\arctan\displaystyle\frac{A_{n}}{B_{n}},$
(16)

where $\delta_{n}=\pm 1$ is selected to make positive $A_{n}\geq 0$ and
$B_{n}>0$.

The very first “remainder” (11) is

$A_{1}=\left|a_{m}\right|,\quad
B_{1}=b_{m},\quad\delta_{1}=\mathop{\hbox{sign}}a_{m}.$ (17)

Let us rephrase (6) separating the fraction in the “remainder” from its sign.
The following equations also give us the expression for $\delta_{n}$ in (16):

$\begin{array}[]{c}\mu_{n+1}\arctan\displaystyle\frac{A_{n+1}}{B_{n+1}}=\arctan\displaystyle\frac{A_{n}}{B_{n}}-\arctan\displaystyle\frac{1}{q_{n}}\\\
\left\\{\begin{array}[]{l}A_{n+1}=\left|q_{n}A_{n}-B_{n}\right|,\\\
\mu_{n+1}=\mathop{\hbox{sign}}(q_{n}A_{n}-B_{n}),\\\
B_{n+1}=q_{n}B_{n}+A_{n},\\\
\delta_{n+1}=\delta_{n}\mu_{n+1}.\end{array}\right.\end{array}$ (18)

Proof. Substitute the “remainder” in (16) according to (18) and write

$\begin{array}[]{ll}\displaystyle\frac{\pi}{4}=m\arctan\displaystyle\frac{1}{q_{0}}+\sum_{k=1}^{n-1}\delta_{k}\arctan\displaystyle\frac{1}{q_{k}}+\delta_{n}\arctan\displaystyle\frac{1}{q_{n}}+\delta_{n}\mu_{n+1}\arctan\displaystyle\frac{A_{n+1}}{B_{n+1}}\\\\[16.0pt]
\phantom{\displaystyle\frac{\pi}{4}}=m\arctan\displaystyle\frac{1}{q_{0}}+\sum_{k=1}^{n}\delta_{k}\arctan\displaystyle\frac{1}{q_{k}}+\delta_{n+1}\arctan\displaystyle\frac{A_{n+1}}{B_{n+1}},\quad\hbox{{\it
Q.e.d.}}\end{array}$

If $A_{n}=0$ in (16) (or (18)), there is no “remainder”, and we have the
identity. Let us assume that $A_{n}>0$. From the Fact 2.2, Corollary 2,
minimum $A_{n+1}$ in (18) happens at

$q_{n}=\left\lfloor\displaystyle\frac{B_{n}}{A_{n}}\right\rceil.$ (19)

Because of the Fact 2.2, $A_{n+1}$ will be

$0\leq A_{n+1}=A_{n}\left|q_{n}-\displaystyle\frac{B_{n}}{A_{n}}\right|\leq
A_{n}/2.$ (20)

Since (20) implies that $A_{n+1}<A_{n}$, the sequence $\\{A_{n}\\}$ is finite.
However, compared to all-positive Machin-like formula from the section 3, its
signed version (10) converges faster, having as few as $O(\ln A_{1})$ terms.
We will present a more meaningful estimation in the next subsection.

###  4.3 Characterization

####  4.3.1 The denominator $q_{1}$ of the second term in (10)

Substitute (17) into (13) and use (3) for double arctangent:

$2\arctan\displaystyle\frac{A_{1}}{B_{2}}=\arctan\displaystyle\frac{2A_{1}/B_{1}}{1-\left(A_{1}/B_{1}\right)^{2}}\leq\arctan\displaystyle\frac{1}{q_{0}}$
(21)

Arctangent is a monotonically increasing function. Taking into account the
trivial corollary of (13)

$\arctan\displaystyle\frac{A_{1}}{B_{2}}\leq\displaystyle\frac{1}{2}\arctan\displaystyle\frac{1}{q_{0}}<\arctan\displaystyle\frac{1}{q_{0}},$

we can write $A_{1}/B_{1}<1/q_{0}<1$. Unequilibrium (21) transforms into

$\displaystyle\frac{2A_{1}}{B_{1}}\leq\displaystyle\frac{1}{q_{0}}\left(1-\displaystyle\frac{A_{1}^{2}}{B_{1}^{2}}\right)<\displaystyle\frac{1}{q_{0}}\quad\Rightarrow\quad\displaystyle\frac{B_{1}}{A_{1}}>2q_{0}.$
(22)

The second term is selected according to (19). Use the Fact 2.2:

$q_{1}=\left\lfloor\displaystyle\frac{B_{1}}{A_{1}}\right\rceil\geq\displaystyle\frac{B_{1}}{A_{1}}-0.5>2q_{0}-0.5,\quad
q_{0},q_{1}\in\mathbb{N}\quad\Rightarrow$ $q_{1}\geq 2q_{0}$ (23)

####  4.3.2 The growth of $\\{q_{n}\\}$

Let us estimate $q_{n+1}$ using its definition (19) and equations (18), (20).
Because of the nearest integer’s properties,

$q_{n+1}\geq\displaystyle\frac{B_{n+1}}{A_{n+1}}-0.5\geq\displaystyle\frac{q_{n}B_{n}+A_{n}}{A_{n}/2}-0.5=2q_{n}\cdot\displaystyle\frac{B_{n}}{A_{n}}+2-0.5.$
$q_{n}\leq\displaystyle\frac{B_{n}}{A_{n}}+0.5\quad\Rightarrow$
$\begin{array}[]{l}q_{n+1}\geq
2q_{n}(q_{n}-0.5)+1.5>2q_{n}^{2}-q_{n}+1\\\\[5.0pt]
\phantom{q_{n+1}}>2q_{n}^{2}-2q_{n}+1=q_{n}^{2}+(q_{n}-1)^{2}\geq
q_{n}^{2}\quad\Rightarrow\end{array}$ $q_{n+1}>q_{n}^{2}.$

Let us recalculate the estimation for $q_{n}$ via $q_{1}$. Lemma.

$q_{n}>q_{1}^{2^{n-1}}\quad\forall\,n\geq 2$ (24)

Proof by induction. At $n=2$, observe:

$q_{n}=q_{2}>q_{1}^{2}=q_{1}^{2^{n-1}}.$

If the Lemma’s statement holds at $n$, verify it at $n+1$:

$q_{n+1}>q_{n}^{2}>\left(q_{1}^{2^{n-1}}\right)^{2}=q_{1}^{2\cdot
2^{n-1}}=q_{1}^{2^{n}},\quad\hbox{{\it Q.e.d.}}$

####  4.3.3 The number of terms $N$ in (10)

Since the overall length $N$ depends on the nominator $A_{1}$, let us express
it first. From (12), notice that

$\left\\{\begin{array}[]{l}a_{k+1}<q_{0}a_{k}\\\
b_{k+1}>q_{0}b_{k}\end{array}\right.\quad\Rightarrow\quad\left\\{\begin{array}[]{l}a_{k}<q_{0}^{k}\\\
b_{k}>q_{0}^{k}\end{array}\right.\quad\Rightarrow\quad b_{k}>a_{k}$

Use (12) again for $b_{k}$:

$b_{k+1}=q_{0}b_{k}+a_{k}<(q_{0}+1)b_{k}\quad\Rightarrow\quad
b_{k}<(q_{0}+1)^{k}$

Taking into account the arctangent bounds, we estimate $m$:

$\begin{array}[]{l}\arctan\displaystyle\frac{1}{q_{0}}=\displaystyle\frac{1}{q_{0}}\cdot\left(1+O\left(\displaystyle\frac{1}{q_{0}^{2}}\right)\right),\quad\hbox{therefore,
from~{}(\ref{remainder-zero}),}\\\\[10.00002pt]
m=\displaystyle\frac{\pi}{4}\left/\arctan\displaystyle\frac{1}{q_{0}}\right.+O(1)=\displaystyle\frac{\pi}{4}\left/\displaystyle\frac{1}{q_{0}}\right.\cdot\left(1+O\left(\displaystyle\frac{1}{q_{0}^{2}}\right)\right)+O(1)\\\\[10.00002pt]
\phantom{m}=\displaystyle\frac{\pi q_{0}}{4}+O(1).\end{array}$

Remember from (17), (22), that $B_{1}=b_{m}$ and $A_{1}/B_{1}<1/(2q_{0})$. Let
us estimate $A_{1}$:

$\begin{array}[]{l}A_{1}=O\left(\displaystyle\frac{(q_{0}+1)^{\pi
q_{0}/4+O(1)}}{2q_{0}}\right)\\\\[10.00002pt]
\phantom{A_{1}}=O\left(q_{0}^{\pi
q_{0}/4}\cdot\displaystyle\frac{q_{0}^{O(1)}}{2q_{0}}\cdot\left(1+\displaystyle\frac{1}{q_{0}}\right)^{\displaystyle\pi
q_{0}/4+O(1/q_{0})}\right)\\\\[10.00002pt] \phantom{A_{1}}=O(q_{0}^{\pi
q_{0}/4})\end{array}$

Overall, taking into account the decay of $\\{A_{n}\\}$ (20), we can estimate
the total length of our Machin-like identity (10):

$N=O(\log_{2}A_{1})=O(\ln q_{0}^{\pi q_{0}/4})=O(q_{0}\ln q_{0}).$ (25)

###  4.4 Lehmer’s measure

Let us summarize what we know so far (23), (24), on the denominators of (10)
and then substitute that into (2).

$n=0:\quad$ $\log_{10}q_{0}$ $n=1:$
$\log_{10}q_{1}\geq\log_{10}2+\log_{10}q_{0}>\log_{10}q_{0}$ $\dots$ $n>1:$
$\log_{10}q_{n}>2^{n-1}\log_{10}q_{1}>2^{n-1}\log_{10}q_{0}$ (26)

Then,

$\begin{array}[]{l}\lambda=\displaystyle\frac{1}{\log_{10}q_{0}}+\displaystyle\sum_{n=1}^{N}\displaystyle\frac{1}{\log_{10}q_{n}}<\displaystyle\frac{1}{\log_{10}q_{0}}+\displaystyle\frac{1}{\log_{10}q_{0}}\cdot\displaystyle\sum_{n=1}^{N}2^{-n+1}\\\\[15.00002pt]
\phantom{\lambda}<\displaystyle\frac{1}{\log_{10}q_{0}}\cdot\left(1+2\displaystyle\sum_{n=1}^{\infty}2^{-n}\right)=\displaystyle\frac{3}{\log_{10}q_{0}}.\end{array}$
(27)

So, the Lehmer’s measure $\lambda$ can be made arbitrarily small by selecting
the appropriate large denominator $q_{0}$ of the first term.

##  5 Partial series

Let us consider the task of calculating $\pi$ with the precision $\varepsilon$
known in advance. The terms in the Machin-like formula (10) decrease very
fast. According to (24) and (25), the last term can be very rougly estimated
as $O(q_{0}^{-q_{0}^{\pi q_{0}/4}})$. The chances are, that all the higher
terms starting with a certain limit will be much less than the required
precision. Let us briefly discuss the strategy for selecting the significant
terms of (10) depending on $\varepsilon$.

For convenience, re-write (10) in form

$\displaystyle\frac{\pi}{4}=mX_{0}+\displaystyle\sum_{n=1}^{N}X_{n},\quad\hbox{where
$X_{n}=\delta_{n}\arctan\displaystyle\frac{1}{q_{n}}$ and
$\left|\delta_{n}\right|=1$.}$ (28)

Let us consider $0<\varepsilon_{1}\ll 1$ such that for some $n$ fraction
$1/q_{n}<\varepsilon_{1}$. Then

$\begin{array}[]{l}\left|X_{n}\right|<1/q_{n}<\varepsilon_{1},\\\
\left|X_{n+1}\right|<\varepsilon_{1}^{2},\\\
\left|X_{n+2}\right|<\varepsilon_{1}^{4}<\varepsilon_{1}^{3},\\\ \dots\\\
\left|X_{n+k}\right|<\varepsilon_{1}^{k+1},\\\ \dots\\\ \end{array}$

Therefore, the difference between the exact identity (28) and the sum of its
beginning terms up to the $X_{n-1}$ will be

$\varepsilon_{\mathrm{tail}}=\left|\displaystyle\sum_{k=n}^{N}X_{k}\right|\leq\displaystyle\sum_{k=n}^{N}\left|X_{k}\right|<\varepsilon_{1}\displaystyle\sum_{k=0}^{\infty}\varepsilon_{1}^{k}=\displaystyle\frac{\varepsilon_{1}}{1-\varepsilon_{1}}.$

To be specific, let us introduce $\varepsilon_{2}$ as function of
$\varepsilon_{1}$, $m$, and $n$:

$\varepsilon_{2}=\displaystyle\frac{\varepsilon_{1}}{(1-\varepsilon_{1})(m+n-1)},$

The significant terms in (28) are computed using the Maclaurin series for
arctangent (9). The error of calculating the arctangent function depends on
the last term $K$ in the arctangent decomposition. Let us require that the
error for each $i$-th term in (28) will be less than $\varepsilon_{2}$. The
length $K_{i}$ of the Maclaurin decomposition shall be large enough to satisfy
the condition:

$\mathrm{error}_{i}\leq\displaystyle\frac{1}{(2K_{i}+3)\,q_{i}^{2K_{i}+3}}<\varepsilon_{2}.$

The actual total error $\tilde{\varepsilon}$ of our computation will be then

$\tilde{\varepsilon}\leq(m+n-1)\,\varepsilon_{2}+\varepsilon_{\mathrm{tail}}=\displaystyle\frac{2\varepsilon_{1}}{1-\varepsilon_{1}}.$

If we select $\varepsilon_{1}$ based on the required precision $\varepsilon$,
we get the bound for the actual error:

$\varepsilon_{1}=\displaystyle\frac{\varepsilon}{2+\varepsilon}\quad\Rightarrow\quad\tilde{\varepsilon}\leq\varepsilon.$
(29)

Notice that all the variables involved in this discussion are functions of
$q_{0}$ and $\varepsilon$, including the length $n$ of the partial Machin-like
formula and the lengths of Maclaurin approximations for each accepted
arctangent term. Selection of $\varepsilon_{1}$ and $\varepsilon_{2}$ as shown
above ensures that the maximum calculation error will be not greater than
$\varepsilon$, so we have constructed the practical method to calculate $\pi$.

Finally, Lehmer’s measure for the incomplete series can be estimated using
inequality (26). The last term of the partial Lehmer’s sum is greater than the
sum of all the trailing Lehmer’s terms. (Assuming that the partial Machin-like
formula has at least $2$ terms.) Thus, doubling the last term in the partial
Lehmer’s sum provides us with the upper bound for the actual Lehmer’s measure.

##  6 Numerical experiment

Computational experiments were conducted for several starting $q_{0}$-s. The
implementation below uses equations (12), (15), (18), and (26). Python 3.6 was
used as the platform with long integer arithmetic support. For partial Machin-
like expressions we set the arbitrary limit of 1 million decimal digits.

    # Computation of recurrent Machin terms. Optionally: partial

    import math;

    def log10(x) :
        return math.log(x) / math.log(10);

    def str_x_or_lg(nm, x) :   # def. lg(x) := log10(x)
        if x < 1e200 :
            return nm + " " + str(x);
        return "lg " + nm + " " + str(log10(x));

    part = int(input("Partial ? (1=yes, 0=no) > "));
    Q = int(input("Start Q > "));

    A = 1;
    B = 1;
    QS = 1;
    M = 0;
    Q0 = Q;
    QQ = [];
    br = 0;

    # first term
    while A*Q >= B :
        print(str_x_or_lg("A", A), str_x_or_lg("B", B));

        M = M + 1;
        tmpA = A * Q - B;
        B = B * Q + A;
        A = tmpA;

    # check other approximation for the first term
    tmpA = A * Q - B;
    tmpB = B * Q + A;
    if A * tmpB + B * tmpA > 0 :
        QS = -1;
        M = M + 1;
        A = -tmpA;
        B = tmpB;
        print(str_x_or_lg("A", A), str_x_or_lg("B", B));

    print("M", M, "\n---");

    # higher terms
    while A > 0 :
        print(str_x_or_lg("A", A), str_x_or_lg("B", B));

        Q = (B + A - 1) // A;

        tmpA = A * Q - B;
        tmpA1 = tmpA - A;

        if abs(tmpA) > abs(tmpA1) :
            Q = Q - 1;
            tmpA = tmpA1;

        B = B * Q + A;
        A = tmpA;

        QQ.append(QS * Q);  # trace the sign in q_k

        if A < 0 :
            A = -A;
            QS = -QS;

        if part > 0 and log10(Q) > 1000000 :
            print("break");
            br = 1;
            break;  # for partial series

    # Lehmer measure
    L = 1 / log10(Q0);
    print("---\nM", M, str_x_or_lg("Q", Q0));
    for q1 in QQ :
        L = L + 1 / log10(abs(q1));
        str_sign = ("(+)" if q1>0 else "(-)");
        print(str_sign, str_x_or_lg("Q", abs(q1)));

    if br > 0 :  # remaining sum in the partial series is less
        L = L + 1 / log10(abs(QQ[-1]));     # than the last term

    # Pi sum for sanity check
    S = M*math.atan(1/Q0);
    for q1 in QQ :
        S = S + math.atan(1/q1);

    say = ("(brk)\n---\nLehm <" if br>0 else "---\nLehm");
    print(say, L, "\nPi", 4*S);

The program takes $q_{0}$ in the input and builds the corresponding Machin-
like identity. For the cases when the running variables become too large,
there is an option to terminate the calculation and make the formula partial.
The program prints the starting multiplier $m$, denominators $q_{k}$, the
signs $\delta_{k}$, and the Lehmer’s measure $\lambda$ in the end. If a
$q_{k}$ is too large for display (we selected the limit of $200$ digits), its
decimal logarithm is printed. For sanity check only, the computed $\pi$ is
also printed.

A few outputs for various $q_{0}$ are listed below.

At $q_{0}=5$, the program prints the original Machin identity:

    M 4 Q 5
    (-) Q 239
    ---
    Lehm 1.851127652316856
    Pi 3.1415926535897936

Let us provide the following identities to demonstrate how to interpret the
Python program output. The numbers can be copy-pasted into a software capable
of scientific computations such as Mathematica to validate the results of our
study.

For $q_{0}=7$:

    M 6 Q 7
    (-) Q 15
    (+) Q 1712
    (-) Q 8886139
    (+) Q 2526830931360443
    ---
    Lehm 2.551666609279759
    Pi 3.1415926535897936

The listing above corresponds to the expression

$\begin{array}[]{l}6\arctan\displaystyle\frac{1}{7}-\arctan\displaystyle\frac{1}{15}+\arctan\displaystyle\frac{1}{1712}\\\\[7.5pt]
\phantom{0}-\arctan\displaystyle\frac{1}{8886139}+\arctan\displaystyle\frac{1}{2526830931360443}\end{array}$

The code to enter to Wolfram Alpha (link: https://www.wolframalpha.com/input)

    Pi/4 == 6 ArcTan[1/7] - ArcTan[1/15] + ArcTan[1/1712]
        - ArcTan[1/8886139] + ArcTan[1/2526830931360443]

For $q_{0}=8$. Notice that while Python supports arbitrarily long integers, it
does not support arbitrary precision floating point numbers:

    M 6 Q 8
    (+) Q 25
    (-) Q 1407
    (+) Q 4150619
    (+) Q 77950325308084
    (+) Q 28355848339635153147414863515
    (-) Q 2412162405181169014685016537064715579879917878585649329193
    ---
    Lehm 2.4159383360928026
    Pi 3.141592653589793

$\begin{array}[]{l}6\arctan\displaystyle\frac{1}{8}+\arctan\displaystyle\frac{1}{25}-\arctan\displaystyle\frac{1}{1407}+\arctan\displaystyle\frac{1}{4150619}\\\\[7.5pt]
+\,\arctan\displaystyle\frac{1}{77950325308084}+\arctan\displaystyle\frac{1}{28355848339635153147414863515}\\\\[7.5pt]
-\,\arctan\displaystyle\frac{1}{2412162405181169014685016537064715579879917878585649329193}\end{array}$

    Pi/4 ==
    6 ArcTan[1/8] + ArcTan[1/25] - ArcTan[1/1407] + ArcTan[1/4150619]
        + ArcTan[1/77950325308084]
        + ArcTan[1/28355848339635153147414863515]
        - ArcTan[1/24121624051811690146850165370
                   64715579879917878585649329193]

For $q_{0}=9$:

    M 7 Q 9
    (+) Q 93
    (+) Q 22055
    (+) Q 5085558009
    (+) Q 767266041127734416424
    (+) Q 1766091533603478722982708121680411788426907
    ---
    Lehm 1.9607629078499424
    Pi 3.1415926535897927

    Pi/4 == 7 ArcTan[1/9] + ArcTan[1/93] + ArcTan[1/22055]
        + ArcTan[1/5085558009] + ArcTan[1/767266041127734416424]
        + ArcTan[1/1766091533603478722982708121680411788426907]

For $q_{0}=10$:

    M 8 Q 10
    (-) Q 84
    (-) Q 21342
    (-) Q 991268848
    (-) Q 193018008592515208050
    (-) Q 197967899896401851763240424238758988350338
    (-) Q 117573868168175352930277752844194126767991
          915008537018836932014293678271636885792397
    ---
    Lehm 1.9473700443296986
    Pi 3.141592653589794

    Pi/4 == 8 ArcTan[1/10] - ArcTan[1/84] - ArcTan[1/21342]
        - ArcTan[1/991268848] - ArcTan[1/193018008592515208050]
        - ArcTan[1/197967899896401851763240424238758988350338]
        - ArcTan[1/117573868168175352930277752844194126767991
                   915008537018836932014293678271636885792397]

Related to the famous approximation of $\pi\approx 22/7$ (see, for example,
[2]), expression for $q_{0}=28$. For the sake of performance, numbers longer
than $200$ decimal digits are represented as their decimal logarithms. While,
strictly speaking, the listing below is no longer an identity, one can
reproduce this test and output exact integer numbers for this Machin-like
formula. Notice that $\log_{10}q_{n}$ almost exactly doubles on each step.

    M 22 Q 28
    (+) Q 56547
    (+) Q 20747394343
    (+) Q 1112172624652580034840
    (-) Q 16659543628852678157467292276729792021493732
    (+) Q 19351587917741573692739048650182250035782554284801229
          80428023197249578624178441690588894
    (+) Q 14718492206740001931852838656976183022784010091410392
          42953147036168205460675285916208006732990521412670908
          69513168086930986444104325857945434713227531064709901
          94861973862674124
    (-) lg Q 350.7305238264204
    (-) lg Q 702.0893561664352
    (+) lg Q 1404.5900031211877
    (+) lg Q 2809.9358190450657
    (-) lg Q 5620.463702225073
    (+) lg Q 11241.25183905937
    (-) lg Q 22484.181013176003
    (-) lg Q 44968.75144493231
    (-) lg Q 89937.82819599868
    (-) lg Q 179876.09422636102
    (+) lg Q 359752.6872249542
    (+) lg Q 719508.3122952792
    (-) lg Q 1439017.5723335177
    (-) lg Q 2878035.9207072803
    (-) lg Q 5756072.228487223
    (-) lg Q 11512146.246898009
    ---
    Lehm 1.091872372535026
    Pi 3.141592653589793

To illustrate how Lehmer’s measure decreases with the first term decrease, the
following is the example with $q_{0}=100000$. Partial mode is on. Notice that
(27) provides a close but somewhat larger bound $\lambda<0.6$.

    M 78540 Q 100000
    (-) Q 544491
    (+) Q 783664894308
    (+) Q 1303088915612811138696591
    (+) Q 7636018810382840305552700218709810164960367081459
    (+) Q 3617852367571966985355152435997587992571879865162469
          49685961215997793793692115231994619217388993130
    (-) Q 1263548033106645763782664751160149068356817307489653
          3051095395918146583179586035763817403459556538268035
          7855957343101646766103100085828143285482786260424809
          4974093639535779278903243063902741055991584
    (+) lg Q 396.72088863680796
    (-) lg Q 793.8533155269043
    (-) lg Q 1588.4525139699301
    (+) lg Q 3177.6648529734907
    (-) lg Q 6356.3408844093965
    (+) lg Q 12713.353728781887
    (-) lg Q 25427.097270768576
    (-) lg Q 50855.26878265154
    (+) lg Q 101710.88560659182
    (+) lg Q 203422.37580891087
    (-) lg Q 406845.1926757058
    (-) lg Q 813691.1891000423
    (-) lg Q 1627383.4447412174
    (brk)
    ---
    Lehm < 0.5405713556036384
    Pi 3.141592653589794

##  7 Supplement

Fact 2.1. Let us provide the formal proof that the single-term identity (1) is
only possible if $q_{0}=1$.

Proof by contradiction. Let $q_{0}>1$, and the single-term identity is
possible. Since $q_{0}>0$, we have $m>0$:

$\arctan\displaystyle\frac{1}{1}-\underbrace{\arctan\displaystyle\frac{1}{q_{0}}-\dots-\arctan\displaystyle\frac{1}{q_{0}}}_{\hbox{$m$
times}}=0.$

The recurrent formula (6) defines the sequences $\\{A_{n}\\}$ and
$\\{B_{n}\\}$ that describe the result of the subtractions. Because of their
definition, $A_{n}$ and $B_{n}$ are polynomials of variable $q_{0}$ with whole
coefficients. By our hypothesis, $q_{0}$ is the root of $A_{m}(q_{0})=0$.

Let us trace how the constant term evolves in $A_{n}$ and $B_{n}$. We have

$\begin{array}[]{l}\left\\{\begin{array}[]{l}A_{n}(q_{0})=a_{n}^{n}q_{0}^{n}+\dots+a_{n}^{0},\\\
B_{n}(q_{0})=b_{n}^{n}q_{0}^{n}+\dots+b_{n}^{0},\end{array}\right.\\\\[10.00002pt]
\left\\{\begin{array}[]{l}A_{0}=1,\\\
B_{0}=1,\end{array}\right.\quad\Rightarrow\quad\left\\{\begin{array}[]{l}a_{0}^{0}=1,\\\
b_{0}^{0}=1,\end{array}\right.\\\\[10.00002pt]
\left\\{\begin{array}[]{l}A_{n+1}(q_{0})=q_{0}A_{n}(q_{0})-B_{n}(q_{0})=a_{n}^{n}q_{0}^{n+1}+\dots-
b_{n}^{0},\\\
B_{n+1}(q_{0})=q_{0}B_{n}(q_{0})+A_{n}(q_{0})=b_{n}^{n}q_{0}^{n+1}+\dots+a_{n}^{0},\end{array}\right.\quad\Rightarrow\\\\[10.00002pt]
\left\\{\begin{array}[]{l}a_{n+1}^{0}=-b_{n}^{0},\\\
b_{n+1}^{0}=a_{n}^{0}.\end{array}\right.\end{array}$

The constant coefficients change according to the table:

$n$ | $a_{n}^{0}$ | $b_{n}^{0}$
---|---|---
0 | 1 | 1
1 | -1 | 1
2 | -1 | -1
3 | 1 | -1
4 | 1 | 1 loops to $n=0$
5 | -1 | 1
$\dots$ | $\dots$ | $\dots$

$a_{n}^{0}$ is either $1$ or $-1$ for any $n$. Since $q_{0}$ is the root of
polynomial $A_{m}$,

$\left(a_{m}^{m}q_{0}^{m-1}+\dots+a_{m}^{1}\right)\cdot q_{0}+a_{m}^{0}=0,$

and the expression in the brackets is a whole number. Therefore, $q_{0}$ is
the divisor of $a_{m}^{0}$. It is impossible if $q_{0}>1$. Q.e.d.

Fact 2.2. $\left|\left\lfloor x\right\rceil-x\right|\leq 1/2$.

Proof. a) If $x\in\mathbb{Z}$, the left part is $0$. Let us consider cases
$x\notin\mathbb{Z}$.

b) $\left|\left\lfloor x\right\rfloor-x\right|<\left|\left\lceil
x\right\rceil-x\right|$ and $\left\lfloor x\right\rceil=\left\lfloor
x\right\rfloor$ (see (8)). Using (7), derive:

$\begin{array}[]{l}x-\left\lfloor x\right\rfloor<\left\lceil
x\right\rceil-x=\left\lfloor x\right\rfloor+1-x\quad\Rightarrow\\\
2\left|x-\left\lfloor x\right\rceil\right|=2x-2\left\lfloor
x\right\rfloor<1.\end{array}$

c) $\left|\left\lfloor x\right\rfloor-x\right|\geq\left|\left\lceil
x\right\rceil-x\right|$ and $\left\lfloor x\right\rceil=\left\lceil
x\right\rceil$. Similar to case (b),

$x-(\left\lceil x\right\rceil-1)\geq\left\lceil
x\right\rceil-x\quad\Rightarrow\quad 1\geq 2\left|\left\lfloor
x\right\rceil-x\right|.$

In all 3 cases, the Fact 2.2 statement holds. Q.e.d.

Fact 2.3 Corollary. $\forall\,x\in(0,1)$, $x-x^{3}/3<\arctan x<x$.

Proof. Let $K=1$ and substitute $\varepsilon$ with its estimation in (9).
Opening the inequality with absolute value,

$-\displaystyle\frac{x^{5}}{5}\leq\arctan
x-\left(x-\displaystyle\frac{x^{3}}{3}\right)\leq\displaystyle\frac{x^{5}}{5},$
$\arctan
x-x\leq-\displaystyle\frac{x^{3}}{3}+\displaystyle\frac{x^{5}}{5}=-\displaystyle\frac{x^{3}}{3}\left(1-\displaystyle\frac{3x^{2}}{5}\right)<0\quad\forall\,x\in(0,1).$

Let $K=2$. Similar to the above,

$\arctan
x-\left(x-\displaystyle\frac{x^{3}}{3}+\displaystyle\frac{x^{5}}{5}\right)\geq-\displaystyle\frac{x^{7}}{7},$
$\arctan
x-\left(x-\displaystyle\frac{x^{3}}{3}\right)\geq\displaystyle\frac{x^{5}}{5}\left(1-\displaystyle\frac{5x^{2}}{7}\right)>0\quad\forall\,x\in(0,1),\quad\hbox{{\it
Q.e.d.}}$

##  8 Conclusion

We presented the method to build the simple Machin-like identity for any
selected starting term. We proved that Lehmer’s measure of our identity can be
made arbitrarily small by selecting a single parameter. We presented the
practical method, and the Python program to compute parameters in our
identity, along with the algorithm for the actual computation of $\pi$ with
pre-defined precision.

##  9 Acknowledgements

I am grateful to Dr. Abrarov, S. M. for the discussion of the results.

##  10 Special thanks

To my wife, Ekaterina, for the encouragement to put my findings on paper.

To the open-source community for the tools to compute and create.

## References

  * [1] S. M. Abrarov, R. Siddiqui, R. K. Jagpal, and B. M. Quine, A new form of the Machin-like formula for $\pi$ by iteration with increasing integers. arXiv:2108.07718v5 [math.GM] 17 Apr 2022.
  * [2] L. Berggren, J. Borwein, P. Borwein, Pi: A Source Book. Springer-Verlag, 1997.
  * [3] D. H. Lehmer, On arccotangent relations for $\pi$. Amer. Math. Monthly 45 (1938), 657-664. DOI: 10.1080/00029890.1938.11990873
  * [4] M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions. Natl. Bureau of Standards, Applied Math. Ser. 55, 10th ed. 1972.
  * [5] R. D. Carmichael, The Theory of Numbers. The Scientific Press, 1914.
  * [6] R. L. Graham, D. E. Knuth, O. Patashnik, Concrete Mathematics: A Foundation for Computer Science. Addison-Wesley, 1994. ISBN 0-201-55802-5
  * [7] M. D. Weir, J. Hass, G. B. Thomas, Thomas’ calculus: early transcendentals. 12th ed. 2009. ISBN 978-0-321-58876-0
  * [8] E. Herman, G. Strang, Calculus Volume 2. Rice Univ. 2017. https://openstax.org/details/books/calculus-volume-2