Astroparticle Physics 62 (2015) 21–29

Contents lists available at ScienceDirect

Astroparticle Physics
journal homepage: www.elsevier.com/locate/astropart

Measurement of the radial density gradient of cosmic ray in the
heliosphere by the GRAPES-3 experiment
H. Kojima a,b, H.M. Antia b,c, S.R. Dugad b,c, S.K. Gupta b,c,⇑, Y. Hayashi b,d, P. Jagadeesan b,c, A. Jain b,c,
S. Kawakami b,d, P.K. Mohanty b,c, T. Nonaka b,d, A. Oshima b,e, B.S. Rao b,c, S. Shibata b,e,
The GRAPES-3 CollaborationaFaculty of Engineering, Aichi Institute of Technology, Toyota City 470-0392, JapanThe GRAPES-3 Experiment, Cosmic Ray Laboratory, Ooty 643 001, India
Tata Institute of Fundamental Research, Mumbai 400 005, IndiadGraduate School of Science, Osaka City University, Osaka 558-8585, JapaneCollege of Engineering, Chubu University, Kasugai, Aichi 487-8501, Japanbca r t i c l e

i n f oArticle history:
Received 24 April 2014
Received in revised form 13 June 2014
Accepted 2 July 2014
Available online 10 July 2014
Keywords:
Cosmic rays
Heliosphere
Solar wind
Convection
Diffusiona b s t r a c t
A radial anisotropy in the ﬂux of cosmic rays in heliosphere was theoretically predicted by Parker andothers within the framework of the diffusion–convection mechanism. The solar wind is responsible forsweeping out the galactic cosmic rays, creating a radial density gradient within the heliosphere. This gradient coupled with the interplanetary magnetic ﬁeld induces a ﬂow of charged particles perpendicular to
the ecliptic plane which was measured and correctly explained by Swinson, and is hereafter referred as‘Swinson ﬂow’. The large area GRAPES-3 tracking muon telescope offers a powerful probe to measure theSwinson ﬂow and the underlying radial density gradient of the galactic cosmic rays at a relatively highrigidity of 100 GV. The GRAPES-3 data collected over a period of six years (2000–2005) were analyzed
and the amplitude of the Swinson ﬂow was estimated to be (0.0644 ± 0.0008)% of cosmic ray ﬂux which
was an 80r effect. The phase of the maximum ﬂow was at a sidereal time of (17.70 ± 0.05) h which was
18 min earlier than the expected value of 18 h. This small 18 min phase difference had a signiﬁcance of
6r indicating the inherent precision of the GRAPES-3 measurement. The radial density gradient of thegalactic cosmic rays at a median rigidity of 77 GV was found to be 0.65% AU1.
Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction
The heliosphere comprises the region of space surrounding the
Sun where the solar wind exerts sufﬁcient pressure to sweep out
the galactic cosmic rays into the interstellar medium. This region
contains the interplanetary magnetic ﬁeld (IMF) which is connected to the Sun, and is also frozen into the solar wind. The
anchoring of the IMF to the solar surface and the ﬂow of the solar
wind forces the IMF outwards into a spiral form due to solar rotation. The solar wind results in the convection of charged particles
away from the Sun, and well beyond the heliosphere lies the interstellar region where the galactic cosmic rays become largely isotropic. The random motion and the interactions of the galactic cosmic
rays cause them to cross the boundary and enter the heliosphere.
⇑ Corresponding author at: Tata Institute of Fundamental Research, Mumbai 400
005, India.
E-mail address: gupta@grapes.tifr.res.in (S.K. Gupta).
http://dx.doi.org/10.1016/j.astropartphys.2014.07.003
0927-6505/Ó 2014 Elsevier B.V. All rights reserved.

These particles gyrate around the IMF but the presence of irregularities in the IMF scatters them away from a regular gyro motion.
Also these IMF irregularities impact the convection of charged particles by mediating the interaction of the solar wind with theincoming cosmic rays. The overall impact of this phenomena is diffusion and convection of the particles from the heliospheric boundary towards the Sun. However, this diffusive motion of particles isinﬂuenced by the gradient and curvature drifts in the IMF [1–3].
The solar modulation in the heliosphere results in streaming ofparticles that produces an anisotropy in the ﬂux of galactic cosmic
rays. This streaming may be decomposed into two major components, ﬁrst one in the ecliptic plane and the second perpendicular
to this plane. The galactic cosmic rays entering the Earth’s atmosphere produce a number of secondary particles. However, among
the secondary particles, neutrons and muons reach the ground leveland may be studied by the detectors placed on the surface of the
Earth. The streaming in the ecliptic plane can be recorded by groundbased detectors as a diurnal variation in the counting rate, super-

22

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

posed on a larger isotropic component. The streaming of primarycosmic rays responsible for this variation is known as the Solar Diurnal Anisotropy [1]. On the other hand, the streaming perpendicular
to the ecliptic plane in the north–south direction is observed as asidereal diurnal variation by the ground based detectors [4].
Several studies using the muon telescopes have shown the presence of such a sidereal variation along the north–south direction. At
energies below 100 GeV, the sidereal diurnal variations caused by
the solar modulation of the galactic cosmic rays in the heliospherebecome dominant. This variation results from the streaming of cosmic rays B  Gr perpendicular to the plane of the ecliptic along thenorth–south direction, here B is the IMF and Gr the radial density
gradient of the cosmic rays in the heliosphere [4]. This streaming
results in generation of a diurnal variation in sidereal time, due tothe tilt of the rotation axis of Earth relative to the plane of the ecliptic. Therefore, a muon telescope observes a maximum counting rate
at 18 h local sidereal time, if the streaming is downward (north–south) and at 6 h local sidereal time if the streaming is upward(south–north). Since Gr points radially outward this streaming process reverses its direction whenever the polarity of the IMF changes.
The region of the heliosphere where the IMF is directed toward theSun, is called the Toward or TW sector, and the streaming is downwards (north–south). Similarly, when the IMF is directed away fromthe Sun, it is termed the Away or AW sector and the streaming is
upwards (south–north). Therefore, the muon telescopes located in
the northern hemisphere of the Earth detect a maximum rate at
6 h in the AW sectors and at 18 h in the TW sectors, in local sidereal
time. This streaming as mentioned earlier is referred as the Swinson
ﬂow. The Swinson ﬂow can be nearly eliminated when averaged
over a year provided the AW and the TW sectors are of similar duration. On the other hand, the characteristics of the Swinson ﬂow may
be extracted by taking the difference of the mean variations in the
AW and TW sectors because the Swinson ﬂow is north–south asymmetric in the ecliptic plane [5].
A spurious sidereal diurnal variation could also be generated byseasonal variation of the solar diurnal anisotropy as explained
below. The line of sight of a muon telescope towards the Sun
undergoes a systematic annual change due to the orbital motionof Earth, resulting in an annual modulation of the solar diurnal variation. This modulation may be resolved into a sidereal and an antisidereal diurnal component of equal amplitude [6,7]. However, thisspurious variation can be easily eliminated by using the anti-sidereal diurnal variation mentioned above [8]. After eliminating thisspurious effect, the surviving effect contains the Swinson ﬂow
and the genuine sidereal diurnal anisotropy. The Swinson ﬂow issensitive to the polarity of the IMF in the vicinity of the Earth, a linear combination of the measurements from the TW and AW sectors may be exploited to extract the Swinson ﬂow from themeasured sidereal diurnal variation as detailed later in this work.
Since the Swinson ﬂow is produced by the solar modulation, it
may be used to obtain the radial density gradient of the charged
cosmic rays [2,9,10].
Various anisotropic phenomena involving the transport of thecosmic rays in the heliosphere, also known as the solar modulation
were theoretically investigated by Parker and others [1]. Thesephenomena are well described by the diffusion–convection model.
In the present work, the data for a relatively long duration of six
years were analyzed. This long duration is expected to result in
the reduction of the effects of the time dependent solar activity.
Thus, a simple diffusion–convection equation without any timedependent terms as expressed by the following equation was used,

Vsw n þ j rn ¼ 0

ð1Þ

where, n and rn, are the density and radial density gradient of cosmic rays, j and Vsw represent the diffusion coefﬁcient and thevelocity of the solar wind, respectively. Here, the radial density

gradient Gr and n are related by Gr ¼ rn=n. The diffusion coefﬁcient
may be used to obtain the value of the mean free path of the cosmic
rays in the heliosphere because these two parameters are two
equivalent descriptions of the same phenomenon. Although, Parker
and others had originally used the diffusion–convection framework
to describe the dynamics of cosmic rays inside the solar system, it
was subsequently employed by several workers to study the acceleration of charged particles in the shock waves produced by supernova explosions ([11] and references contained therein). This
clearly indicated that the experimental studies using the frameworkof diffusion–convection model could be an effective tool to probe
the mechanisms of the cosmic ray acceleration in our galaxy besides
the solar system.
A number of early measurements were made to experimentally
determine different parameters in the Parker equation. However,most of these measurements were made using relatively low
energy cosmic rays observed by the ground based neutron monitors at energies [10 GeV or by the satellite based detectors at evenlower energies [10]. The GRAPES-3 experiment [14] contains alarge area (560 m2) tracking muon telescope [15] that allowed
even a small variation in the muon ﬂux to be measured with highstatistical precision. In the present work, an attempt was made toexperimentally measure the radial density gradient of the galactic
cosmic rays, without recourse to Monte Carlo simulations, and byusing higher energy (77 GeV) galactic cosmic rays. Here, it needs
to be emphasized that only the GRAPES-3 muon data were usedin the present study for estimating the density gradient of galactic
cosmic rays near the Earth. The use of data from an experiment
located at a single geographic location avoided the systematic
uncertainties that are normally associated with the use of datafrom multiple locations due to varying physical conditions and different characteristics of detectors used.

2. GRAPES-3 experimental setup
The GRAPES-3 experiment was set up at Ooty (11.4°N latitude,
76.7°E longitude, and 2200 m altitude) in south India. The GRAPES3 experiment consists of two major components, the ﬁrst component is an extensive air shower array of compact conﬁguration
with an inter-detector separation of 8 m between the adjacent
plastic scintillator detectors, each 1 m2 in area [12,13]. These
detectors were deployed in a symmetric hexagonal geometry. At
present, the array is being operated with about 400 detectors.
The array was designed to study the primary cosmic rays from
10 TeV to 100 PeV and to make precision measurement of theirenergy spectrum and composition [14,16].
The second major component of the GRAPES-3 experiment is alarge area tracking muon telescope that is used to measure themuon content of the cosmic ray showers [15]. This information is
critically important in determining the composition of the primarycosmic rays [17]. This unique instrument was also used to measure
the variation of the cosmic radiation due to the solar activity [18–20]. The muon telescope is capable of providing a high statistics,
directional measurement of the muons. The GRAPES-3 muon telescope consists of a total of 16 individual modules, each of themwith an area of 35 m2. The basic element of the muon telescopeis the proportional counter (PRC) fabricated from a rugged
600 cm long seamless steel pipe with 2.3 mm thick wall and asquare cross-section of 10  10 cm2. A module of sensitive area
35 m2 consists of 232 PRCs arranged in four layers of 58 PRCs each,with alternate layers aligned in mutually orthogonal directions asshown schematically in Fig. 1. Two successive layers of the PRCs
are separated by a 15 cm thick reinforced cement concrete layer
assembled from blocks, each 60 cm  60 cm  15 cm in size. The
four-layer PRC conﬁguration of the muon module allowed

23

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

Layer−3
Layer−2
Layer−1
Layer−0Fig. 1. A schematic of the 4-layer tracking muon telescope module with 58 PRCs per layer. Four layers of PRCs labeled Layer-0, Layer-1, . . . embedded in concrete. Inclined
lines indicate ﬁve parallel muons.

reconstruction of each muon track in two mutually orthogonalplanes. The vertical separation of the two PRC layers in the same
projection plane was 50 cm allowing measurement of the muontrack direction with an accuracy of [6° in each projected plane.
The four layers of PRCs, namely, layers-0, 1, 2, and 3 in mutuallyorthogonal directions were interleaved with absorber concrete
blocks that are clearly visible Fig. 1. Also shown in Fig. 1 are ﬁve
tracks of inclined parallel muons.
To achieve an energy threshold of 1 GeV for vertical muons, anabsorber of thickness 550 g cm2 in the form of concrete blocks
was used as shown schematically in Fig. 1. The concrete blocks
were arranged in the form of an inverted pyramid to shield the
PRCs with absorber coverage up to 45° for incident muons. Consequently, the muon telescope had an energy threshold of sec (h) GeV
for the muons incident at a zenith angle of h. As shown in Fig. 2, thearrival direction of a muon was determined for each triggered PRC
in the lower layer and binned into 13 directions based on the speciﬁc location of the PRC triggered in the upper layer from amongthe 13 PRCs, one directly above (central PRC) and six each on eitherside of the central PRC. This directional binning was carried out ineach of the two orthogonal projection planes, thereby generating a
map of 13  13 = 169 solid-angle directions of the muons as shownin Fig. 3. The contents of the 169 directions were recorded once
every 10 s, thereby generating a continuous record of the directional ﬂux of muons in the sky.
However, to obtain a higher statistical precision, the 169 directions were combined into nine coarse directions as shown schematically in Fig. 3. This was done by combining either a set of3  5 or 5  5 solid-angle directions with the exception of the vertical direction, where the central 3  3 solid-angle directions were

−1

0

1

2

4

6

6
5
4

NW

N

NE

W

V

E

SW

S

SE

−4 −3 −2 −1

0

3
2
1
0
−1
−2
−3
−4
−5
−6
−6 −5

1

2

3

4

5

6

Fig. 3. A schematic of 169 muon directions that were subsequently combined intonine coarse directions. First, 3  3 vertical direction V, four 3  5 central (N, E, W, S),
and four 5  5 outer directions (NE, SE, SW, NW).

combined. This speciﬁc combination scheme resulted in a relatively similar solid-angle coverage for these nine coarse directions.
This choice of the angular segmentation also reduced the dissimilarity in the number of muons recorded in each coarse direction,
because of the comparatively larger ﬂux of muons for the near central directions (N, E, W, S) when compared to the outer directions(NE, SE, NW, SW). The cutoff rigidity at Ooty was 17 GV in the vertical direction and it varied from 14 to 32 GV across the ﬁeld ofview (FOV) of the muon telescope. However, in the present work
only the data from the E, V, and W directions were utilized.
3. Data analysis
3.1. Data cuts

50ο

39 ο

22ο 11ο 0ο −11οFig. 2. Schematic view of muon arrival angle selection based on triggered PRCs, one
in the lower and one from among the 13 PRCs in the upper layer. Triggered PRCs
shown as ﬁlled squares.

As described earlier, the one-dimensional Parker convection–diffusion framework provides a relationship among the radial density gradient of the cosmic rays, the solar wind velocity (Vsw), andthe cosmic ray diffusion coefﬁcient (j). The GRAPES-3 experimentdetected the muons of energy above 1 GeV, that correspond to a
median energy of 66 GeV for the primary cosmic ray protons in
the vertical direction and varied from a low of 64 GeV–92 GeV
across the FOV of the GRAPES-3 muon telescope. The medianH. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

energy of all primary cosmic ray protons within the FOV of muon
telescope was 77 GeV. These values of the median proton energies were estimated with the help of the international geomagnetic
reference ﬁeld model (IGRF-11) of the geomagnetic ﬁeld [21] andMonte Carlo simulations using the CORSIKA package [22]. Thus
the measurements using the GRAPES-3 data described below correspond to a median rigidity of 77 GV. Since the GRAPES-3 tracking muon telescope recorded about 4  109 muons every day or1.5  1012 muons per year, it thus provided a directional surveyof the sky with exceptional statistical precision.
The solar diurnal variation at a rigidity of 77 GV was studied
from the measured ﬂux of muons by adopting the following procedure. The muon data for intervals of one hour were individuallysegregated for each of the nine directions after correcting for thedead-time. Next, a two step process was employed to effectively
identify and remove the gaps and fast transients present in thedata. This was necessary because such transients and data gaps
had the potential to disturb the measurement of the Swinson ﬂow
that varies slowly with time. The root mean square (rms) deviation
labeled r1 for each of the nine directions were calculated. For eachdirection the hourly data set corresponding to a deviation  10r1
were rejected as abnormal. This procedure was repeated a second
time for the data after the ﬁrst step to recalculate the mean and therms (r2 ). Next, the data with a deviation of  5r2 were also
rejected. The implementation of this rigorous two step cut eliminated poor quality data segments that included data gaps, gain
variations, and the transient phenomena such as the Forbush
decrease events of magnitude larger than 3%. Thereafter, the hourly
muon rate for each direction was corrected by using the pressure
coefﬁcients for the corresponding direction for variation in theatmospheric pressure. Finally, the hourly muon rates wereexpressed as fractional changes (in %) from the mean of the data
for the six year period (2000–05) for each of the nine directions.
After completion of the above quality cuts, short-term transients were also removed from the data while retaining the longterm variations including the 27-day, the seasonal (annual) andthe 11-year solar effects and other short-term periodic behaviors
such as the sidereal and solar diurnal variations. Next, to probe
the short-term periodic variations, it became essential to ﬁlter
out the slower long-term variations listed above. The effects oflong-term variation were removed by devising an indirect highpass ﬁlter. The high-pass ﬁlter was implemented by the following
procedure. A running mean H(low)i of the hourly muon rates hi wascalculated for each hour by taking the data of 12 h on either side ofa given hour for each of the nine directions. This running mean
provided an effective low-pass ﬁlter H(low)i of the hourly rates.
Then, the running mean H(low)i was subtracted from the corresponding observed hourly rate hi to obtain the hourly high-passdata H(high)i for the nine directions. For further analysis of thesolar and sidereal diurnal variations, we selected only the days that
contained data for each of the 24 h without any gap.
In order to remove the residual effects of the transient phenomena such as a Forbush decrease event, the calender days with signiﬁcant transient activity were identiﬁed by using the Kiel neutron
monitor data [23]. Due to a lower cutoff rigidity, the Kiel neutron
monitor recorded a large amplitude even for Forbush decrease
events of relatively small magnitude and thus could be used to
identify such events. The following criteria were used to identify
the days of transient activity, (a) daily rate N i for the ith day wascompared with the mean daily rate N M for the previous three days.
Here, the rates N i and N M were expressed in %. If the value ofðN M  N i Þ ¼ DN exceeded 4% then the data for the ith day for allnine directions were attributed to be inﬂuenced by a transient phenomena and rejected, (b) The days following the ith day deﬁned
above were also attributed to the recovery phase of the transientphenomena. If the deviation from N M of the neutron monitor rates

N iþ1 ; N iþ2 ; . . . for the following days exceeded 0:3  DN, then thedata for ði þ 1Þth, ði þ 2Þth, . . . were deemed to lie in the recovery
phase and rejected.
3.2. East–West method
The GRAPES-3 muon data were inﬂuenced by several effects
including those induced by the atmospheric pressure and temperature. The variation in muon rates due to pressure changes werecorrected as described earlier. However, it was not possible to correct for the effects of temperature and other unknown factors.
Therefore, we adopted a powerful technique proposed by Kolhörster and others [24–26], which is capable of eliminating the effects
listed above. This East–West method was implemented as
described below.
The variation in the intrinsic anisotropy in space may be calculated by using the measurements from the E and W directions ands, where s is the time offset between E and W directions as shown
below,

DI ¼ Iðt þ DtÞ  IðtÞ ¼ Iðt þ s=2Þ  Iðt  s=2Þ ¼ NE ðtÞ  NW ðtÞð2Þ

If the intrinsic anisotropy IðtÞ in space is superposed on an isotropic
component B and a common time dependent variation CðtÞ of atmospheric origin, then the time variation recorded along the V, E, and
the W directions, N V ðtÞ, N E ðtÞ, and N W ðtÞ, respectively, may be written as,

NV ðtÞ ¼ IðtÞ þ B þ CðtÞ

ð3ÞNE ðtÞ ¼ Iðt þ s=2Þ þ B þ CðtÞ

ð4Þ

NW ðtÞ ¼ Iðt  s=2Þ þ B þ CðtÞ

ð5Þ

NE ðtÞ  NW ðtÞ ¼ Iðt þ s=2Þ  Iðt  s=2Þð6Þ

By summing the difference N E ðtÞ  NW ðtÞ obtained from the
GRAPES-3 observations, the time independent anisotropy in space
was reconstructed using the East–West (E–W) method as detailed
below. However, one needed to experimentally measure the timeoffset s corresponding to the hour angle difference between the E
and W directions. The hourly muon rates expressed as percent deviation from the six year mean for E, and W directions for all six yearsof data were each folded modulo 24 h. The time offset s due to thehour angle difference between E and W was precisely calculated byusing the cross correlation of the muon rates from these two directions. The result of this correlation study is shown in Fig. 4. The
maximum correlation between E and W was obtained for an offset

100
90
Correlation (%)

24

80
70
60
50

-12

-9

-6

-3

0

3

6

9

12

Offset (hours)
Fig. 4. Correlation between E and W as a function of time offset (h). Maximum
correlation of 89% observed for offset of 3 h consistent with the hour angle of 3 h
between E and W directions.

25

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

of s = 3 h, exactly as expected from the geometry of the GRAPES-3
muon telescope.
The hourly muon rates expressed as percent deviation from thesix year mean for E, W, and V, directions after folding modulo 24 hfor six year interval are shown in Fig. 5(a)–(c), respectively. The
plot shown in Fig. 5(d) was obtained by using the E–W method
as described below. This was achieved by using the offset of 3 hobtained earlier and the hourly rate for the E and W directions.
Thereafter, the expected hourly time variation for the V direction
was obtained as shown in Fig. 5(d).
3.3. Data classiﬁcation by IMF
Next, we describe the methodology followed in identifying thedays that were characterized by the IMF oriented either toward(labeled ‘TW’) or away (‘AW’) from the Sun, respectively. The geocentric solar ecliptic (GSE) coordinate system with the Earth at theorigin was used. The positive X-axis was deﬁned along the line
joining the Earth to the Sun, and the negative Y-axis was tangent
to the orbit of the Earth along the direction of its motion. The

0.5

(a)

0.0

0.5

(b)

Amplitude (%)

0.0

0.5

(c)

0.0

Table 1Category of IMF directions.

0.5

(d)

0.0

0

positive Z-axis was along the north ecliptic pole forming a righthanded Cartesian coordinate system. The magnetic ﬁeld data
showed that the magnitude of the ‘z’ component of IMF, Bz wasmuch smaller than the other two components, namely, Bx and By
and hence only the latter two components were used [27]. A running mean of the hourly values of Bx and By were calculated over aninterval of 72 h centered on the current hour and tabulated into
four categories as listed in Table 1. As mentioned earlier, the median energy of the primary protons responsible for the muonsdetected in the GRAPES-3 muon telescope was 77 GeV. The
choice of 72 h was dictated by the fact that the Larmor radius ofa 77 GeV proton is 0.36 AU in an observed mean IMF of 4.7 nT.
For a typical velocity of 400 km/s, the solar wind takes about
36 h to travel a distance of 0.36 AU. Thus a 72 h interval centeredon the hour in question allowed a reasonable estimate of the effective IMF to be obtained. Using the above procedure the runningmean of the IMF components were calculated and based on thesign of Bx and By , the data were assigned to one of the four categories listed in Table 1 for each of the 24 h in a calender day.
If the IMF data for each hour in a day fell exclusively into either
the AW or the TW category, then that day was selected for further
analysis. Based on the complete days collected during 2000–05,
about 79% of days fell into either AW or TW categories, and theremaining 21% into the mixed category and hence were excluded
from further analysis as mentioned in Table 2.
It was difﬁcult to directly extract the sidereal diurnal variation
because its amplitude was about an order of magnitude smaller
than the solar diurnal variation. Therefore, we indirectly estimatedthe amplitude of the sidereal diurnal variation as explained below.
Since the sidereal and the solar years synchronize once in a year,
due care was taken to ensure that the data for the whole years
were used to extract these anisotropies. If the solar diurnal anisotropy had a constant amplitude, then the data analyzed for the sidereal diurnal anisotropy for a complete year would wipe out thesolar diurnal anisotropy. However, the small inclination angle of
7.3° between the solar equatorial plane and the plane of the ecliptic introduced a seasonal variation which manifested as a sidereal
and an anti-sidereal variation of equal amplitude. In addition, there
might have been other causes of seasonal variation that could
mimic a sidereal and an anti-sidereal variation. Therefore, an analysis strategy was devised to identify and largely eliminate theseseasonal effects. Since there were gaps in the data due to rejectionof the days when the IMF did not fall into AW or TW category, it
was important to use a method that would not have introduced
any side bands.

#

Bx

By

Type

Action

1
2
3
4

<0
60
P0
>0

<0
P0
60
>0

Mixed
Away
TowardMixed

RejectAcceptAcceptRejectTable 2
Data categories.

6

12
18
IST (Hours)

24

Fig. 5. Percent deviation from hourly mean muon rate for 2000–05 data for, (a) E,(b) W, (c) V, and (d) VEW.

Class

Days

Total days (2000–05)
Full days (All)
Full days (TW)
Full days (AW)Full days (Mix)

2192
2015
784
810
421

26

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

The data for a given calender month, say January 2001 were
sorted into two groups labeled, Jan01AW and Jan01TW containing
the data corresponding to AW and TW categories, respectively.
Next, the mean of the data for each hour for all days present in
Jan01AW were calculated to extract the variation in 1 h intervals
for this category. The same exercise was repeated for the Jan01TW
data. These monthly data pairs were then used for further analysis.
The monthly data for six years (72 sets each) were individuallycombined to yield the TW and AW proﬁles. The mean of the TW
and AW proﬁles largely caused by the solar diurnal anisotropy
for the period 2000–05 are shown in Fig. 6(a) and in Fig. 6(b),respectively. The height of the peaks for the TW and AW categories
as seen from Fig. 6(a) and (b) are very similar clearly indicatingthat the polarity of the IMF had very little effect on the amplitudeof the solar diurnal anisotropy.
In principle, the TW and the AW data might have contained thecontributions from the IMF polarity independent effects such as
the solar diurnal anisotropy and possible IMF polarity dependent
effects. The sum (TW+AW)/2 cancels out the IMF polarity dependent effects, while retaining the IMF polarity independent solar
diurnal anisotropy. The mean solar diurnal variation representedby (TW+AW)/2 for the period 2000–05 when calculated was nearly
identical to the plots shown in Fig. 6(a) and (b). Due to their similar
proﬁles the data for the AW and the TW categories were combined
and the amplitude of the solar diurnal variation was estimated tobe 0.3% which peaked at 13 h Indian Standard Time (IST) [28].
This corresponded to an asymptotic direction of 16 h IST. Similarly,
the quantity (TWAW)/2 canceled out the IMF polarity independent solar diurnal anisotropy, while retaining the IMF polarity
dependent effects. However, this exercise did not show a clear
indication of an IMF dependent anisotropy in the GRAPES-3 data.
A more sensitive analysis is needed to probe the potential presenceof such an IMF polarity dependent solar anisotropy.

By employing the E–W method outlined earlier, the variation in
the TW and the AW sectors were obtained by folding the monthly
data modulo the solar day. Next, the quantity (TWAW)/2, which
is IMF polarity sensitive was calculated for each month. To estimate the sidereal anisotropy the (TWAW)/2 data from successivemonths were shifted backwards by 2 h relative to the previous
month, to synchronize the data to the sidereal clock. The variationof (TWAW)/2 for the period 2000–05 is shown in Fig. 7(a) which
represents the effect of the IMF polarity dependent sidereal anisotropy as measured using the GRAPES-3 muon data.
The presence of the seasonal effects on the solar diurnal variation through an interference introduces both a sidereal and ananti-sidereal variation of equal amplitude. The procedure outlinedin the previous paragraph when repeated by shifting the data forward by 2 h synchronized the data to an anti-sidereal clock, thus
canceling the sidereal variations while retaining the contributionof the seasonal effects as shown in Fig. 7(b). In view of a relatively
small amplitude the contribution of the seasonal effects on theobserved sidereal variation is also expected to be very small. Next,
the anti-sidereal contribution from Fig. 7(b) was used to correct
the measured sidereal variation shown in Fig. 7(a) by using theprescription of Nagashima and coworkers [8] to obtain the correct
amplitude of sidereal variation. This variation representing theSwinson ﬂow is shown in Fig. 7(c). The fact that the amplitude of
the Swinson ﬂow before and after the correction was nearly sameis due to a very small amplitude of the seasonal correction
displayed in Fig. 7(b).

0.1

(a)

0.0

0.5

(a)
-0.1

Amplitude (%)

0.1

Amplitude (%)

0.0

(b)

0.0

0.5
-0.1

(b)

0.1

(c)

0.0

0.0

-0.1

0

6

12

18

24

IST (Hours)
Fig. 6. Muon rate variation (%) due to solar diurnal anisotropy (2000–05) for, (a)TW, (b) AW categories.

0

6

12

18

24

Local SiderealTime (Hours)
Fig. 7. Muon rate variation (%) due to sidereal diurnal anisotropy for 2000–05, (a)(TWAW)/2, (b) anti-sidereal contribution (c) corrected (TWAW)/2

27

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

The variation due to the Swinson ﬂow represented by the corrected (TWAW)/2 shown on Fig. 7(c) displays a minimum at
6 h sidereal time. This may be understood with the help of a schematic of the interplanetary medium that includes the Sun and theEarth as shown in Fig. 8. The gradient in the radial density of thecosmic rays results in radially inward diffusion of cosmic rays.
When the IMF points toward the Sun as shown in Fig. 8(a), the Lorentz force exerted on the cosmic rays diffusing inward induces anorth–south ﬂow, perpendicular to the ecliptic plane that displays
a minimum at 6 h sidereal time and a maximum at 18 h as seenfrom Fig. 8(a). However, when the IMF points away from the Sunthen the resultant ﬂow displays a maximum at 6 h sidereal timeand a minimum at 18 h as shown in Fig. 8(b). Therefore, theparameter (TWAW)/2 provides a more precise measure of thisnorth–south TW ﬂow by eliminating the common sources of errors
caused other phenomena. The symmetry of the plot in Fig. 7(c)implies that the Swinson ﬂow along both, the north to south andthe south to north directions had nearly the same magnitude. Since
the variation shown in Fig. 7(c) displayed a periodic behavior it
was decided to perform a ﬁt by using sine and cosine terms as follows. The data for each hour as ti , and corresponding amplitude
and error as N i , and ri , where i ¼ 1; 2; 3; . . . n ðn ¼ 24Þ. The following function Fðt i Þ was ﬁtted to the data by least squares method,





ti
ti
þ b sin 2pFðti Þ ¼ a cos 2p
24
24

ð7Þ

where n = 24 is the number of data points and the values a and b
were determined as described below. The weighted sum (S) of thesquares of residuals (Di ) is expressed as,

S¼

n
X
D2
i¼1

,

i
2
i

r

n
X
1
i¼1

ð8Þ

r2i

(a)

Toward

B

Sun

18h

Ecliptic plane

0h 6h
18h

Gr

Away

Sun

18h

(b)

Ecliptic plane

6h

S¼




2
n
n 
X
X
ti
ti
 b sin 2p
D2i ¼
Ni  a cos 2p
24
24
i¼1
i¼1

ð9Þ

By minimizing the quantity S with respect to a and b, one obtains,a ¼ 0:0050  0:0052%

ð10Þ

b ¼ 0:0642  0:0052%

ð11Þ

and this ﬁt yielded v ¼ 0:468 for a total of ðn  2Þ ¼ 22 degrees of
freedom which is too small a value to be realistic. Such a small valueof v2 clearly indicated the presence of a large systematic variation,and therefore, the true error was estimated from the sum of squareof the deviation S of the data points from the function listed in Eq.
(7). This was implemented by minimizing the rms deviation of thedata points from the function described in Eq. (7). A proper propagation of the errors led to the following estimate of the rms deviation r,
2

vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
u
n
u 1 X
r¼t
D2
n  2 i¼1 i

ð12Þ

solution of the above equation yields the value of r to be,sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
0:00014694% ¼ 0:0026%
r¼
ð24  2Þ

ð13ÞThis value of r = 0.0026% is approximately seven times smaller thanthe value of 0.018% estimated from the rms deviation of hourly values for 72 months of data. Clearly, the Swinson ﬂow despite showing a sizable variation over a time scale of six years had aremarkable underlying stability that had manifested in a rather stable proﬁle characterized by a small value of r ¼ 0:0026%. By incorporating this new value of r the errors in amplitudes a and b were
estimated to be,a ¼ 0:0050  0:0008%

ð14Þ

b ¼ 0:0642  0:0008%

ð15Þ

From these two parameters the actual amplitude A and the phase /was found to be,

6hB

Since in the present case variation in the value of ri was negligibly
small, the above calculation was carried out with a constant value ofri ¼ r for all i as given below,

18h
0h 6h

Gr

Fig. 8. Schematic of Swinson ﬂow at Earth indicated by colored arrows due todiffusion of cosmic rays in IMF. Curved thick red arrow indicates excess ﬂow andblue thin arrow deﬁcit ﬂow. IMF relative to Sun in, (a) TW sector, sinusoidal ﬂow
with maximum at 18 h, (b) Away sector, sinusoidal ﬂow with maximum at 6 h.
(TWAW)/2 removed common systematic effects and yielded actual TW ﬂow.

A ¼ 0:0644  0:0008%

ð16Þ

w ¼ 17:70  0:05 hð17Þ

Thus the amplitude of the Swinson ﬂow as measured by GRAPES-3
as shown in Fig. 7(c) was (0.0644 ± 0.0008)% with a phase of(17.70 ± 0.05) h, which was close, but not identical to its expectedvalue of 18 h. The difference in the expected and measured phasesof (18 ± 3) min had a signiﬁcance of nearly 6r.
In the present work the direction from south to north, perpendicular to the ecliptic plane was referred as the upward and theone from north to south as the downward direction. The magnitude of the Swinson ﬂow was estimated as explained in the following. It was assumed that the Swinson ﬂow caused an anisotropy ofmagnitude n%, perpendicular to the ecliptic plane along the downward direction. This angular dependence may be expressed as
1 þ n cos h, where h is the angle relative to the perpendicular
described above. A muon telescope on Earth would observe thisanisotropy as a sidereal diurnal variation and for a location on
the equator an amplitude of n sin / would be observed with a maximum at a sidereal time of 18 h, where / is the angle between theequatorial and ecliptic planes. However, a telescope located at alatitude k in the northern hemisphere would observe an amplitudeof n sin / cos k with a maximum also at a sidereal time of 18 h.

28

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

The angle between the upward perpendicular to the ecliptic andthe vertical direction at a latitude k on Earth may be expressed ash1 ¼ p=2  k  / at a sidereal time of 18 h and h2 ¼ p=2  k þ / at
6 h. Therefore, a telescope located at a latitude k in the northernhemisphere observing the cosmic rays in the vertical direction
would observe the difference between the maximum and minimum of the anisotropy, nðcos h1  cos h2 Þ ¼ nfsinðk þ /Þ  sinðk  /Þg ¼ 2n cos k sin /. Thus the magnitude of the anisotropy n(%) may be calculated from the amplitude A (%) estimated earlier
directly from data,n¼

A
cos k sin /ð18Þ

By substituting the values k and / in the above equation,

A ¼ 0:0644  0:0008%


k ¼ 11:4/ ¼ 23:5



ð19Þ
ð20Þ
ð21Þ

The magnitude of the anisotropy caused by Swinson ﬂow was
obtained,

n ¼ 0:1648  0:0020%

ð22Þ

4. Discussion
The large area GRAPES-3 tracking muon detector allowed measurement of the ﬂux of >1 GeV muons along nine directions labeled
NE, N, NW, E, V, W, SE, S, SW, each with a solid angle coverage of0.23 sr. The muons were detected at a rate of 1.8  108 h1enabling small variations in its ﬂux to be reliably measured. A
two step quality cut was imposed on the muon data to remove
the data that had gaps, gain variation of the proportional counters,
and was inﬂuenced by transient phenomena such as the Forbush
decrease events etc. Despite these stringent cuts only 177 d (8%)
out of 2192 d during the six years (2000–05) were rejected. By
the use of the powerful E–W method, we could largely eliminate
contributions due to various systematic effects including the onesof atmospheric origin, and reconstruct the actual time variation
along the V direction. This variation represented the solar diurnal
anisotropy that peaked at 13 h IST locally and 16 h IST
asymptotically.
Using the IMF components Bx, By, Bz [27], the orientation of the
IMF with respect to Sun was determined every hour. Depending on
the orientation of the IMF being either toward (TW) or away (AW)from the Sun during a given day, the data for that day were accordingly classiﬁed as belonging to the TW or AW category. However,during the 24 h in a day when the orientation of the IMF was not
exclusively TW or AW, then that day was classiﬁed as a mixed
day and the data for such days were rejected. A total of 421 out
of 2015 d fell into the mixed category, thereby resulting in a further rejection of about 21% of the data as listed in Table 2. Finally,
784 (TW) and 810 (AW) for a total of 1594 d out of 2015 d survived
the ﬁnal cuts. This data set has formed the base for the analysis
presented here.
The difﬁculties in directly extracting the sidereal diurnal variation in the GRAPES-3 data in the presence of various atmosphericeffects and a much larger solar diurnal variation were overcome
by eliminating these effects using the E–W technique. Since thesidereal and solar years synchronize once in a year the data for
integral number of years were used to extract the sidereal anisotropy. The equal amplitude sidereal and anti-sidereal variations
produced by various seasonal effects were used to eliminate thoseseasonal effects. The data were carefully examined to ensure thatthe gaps in the data due to rejection of the days when the IMF

was in the mixed category did not result in the introduction of
new periodic variations.
The daily data from each calender month were subdivided into
two groups corresponding to the AW and TW categories. Thereafter, the hourly rate for the days in a given month were used to calculate the variation for these two categories to generate monthly
data for six years (72 sets each) for the TW and AW categories.
The six year mean TW and AW proﬁles due to the solar diurnal
anisotropy for the period 2000–05 displayed in Fig. 6(a) and (b)showed that the IMF polarity did not have a large effect on theamplitude of the solar diurnal anisotropy. The mean amplitude of
the solar diurnal anisotropy for 2000–05 was estimated to be0.3%. The anisotropy as shown in Fig. 5(d) peaked at 13 h IST consistent with our measurements for 2006 using a different method[28].
Turning to the main objective of the present work, we next
focused on the sidereal variation. The IMF polarity sensitive quantity (TWAW)/2 was calculated for each month. For obtaining thesidereal variation, the (TWAW)/2 data from successive months
were shifted backwards by 2 h relative to the previous month,
to synchronize the data to the sidereal clock and averaged. The
variation of (TWAW)/2 for the period 2000–05 thus obtained
is shown in Fig. 7(a). The presence of the seasonal effects on
the solar diurnal variation through interference introduces botha sidereal and an anti-sidereal variation of equal amplitude. The
procedure outlined above was repeated for (TWAW)/2 by shifting the data forward by 2 h to synchronize to the anti-sidereal
clock. The anti-sidereal variation thus obtained represented thetiny amplitude of the seasonal effects as shown in Fig. 7(b). Thisvariation was used to correct the observed sidereal variation. The
corrected sidereal variation of (TWAW)/2 is shown in Fig. 7(c),which was the IMF polarity dependent sidereal diurnal variation
representing the projection of the Swinson ﬂow along thedirection V in the FOV of the GRAPES-3. The amplitude of
the sidereal variation before and after the correction was nearly
the same due to a very small amplitude of the correction asshown in Fig. 7(b).
The IMF polarity dependent sidereal diurnal variation shown on
Fig. 7(c) displayed a minimum at 6 h sidereal time, indicating thedirection of the ﬂow to be from south to north, perpendicular to
the ecliptic plane. Similarly, the maximum at 18 h sidereal time
was indicative of the ﬂow from north to south. Since the angle
between the Swinson ﬂow along the north–south (N–S) and the
direction V in GRAPES-3 FOV is minimum at 18 h, the IMF polarity
dependent sidereal diurnal variation projected along V was maximum. Similarly, the angle between S–N Swinson ﬂow and V wasminimum at 6 h, and the amplitude of the Swinson ﬂow projectedalong V was maximum. The symmetry of the plot in Fig. 7(c)implied that the IMF polarity dependent sidereal diurnal variation
along both the N–S and S–N directions was of nearly the same
amplitude. Based on these measurements the amplitude of theSwinson ﬂow was estimated to be (0.0644 ± 0.0008)%.
As described earlier, if the relative magnitude of the Swinson
ﬂow at a rigidity P was n%, the direction of this ﬂow would be perpendicular to the ecliptic plane. The maximum intensity wasobserved if the muon telescope was pointing either along
upstream or downstream of the Swinson ﬂow. The rotation axisof the Earth is inclined by / = 23.5° relative to the ecliptic andthe latitude of the GRAPES-3 site is k = 11.4°. Since the GRAPES-3
is located in the northern hemisphere, thus the maximum southward Swinson ﬂow of n cosðk þ /Þ was expected to be observed
at 18 h, and a minimum of n cosðk  /Þ at 6 h sidereal time. Thus
the amplitude of the IMF polarity dependent sidereal variation
observed along V by the GRAPES-3 telescope is given bynðcosðk þ /Þ  cosðk  /ÞÞ=2 which equals n cos k sin /. Using theobserved value of n cos k sin / of (0.0644 ± 0.0008)% the relative

H. Kojima et al. / Astroparticle Physics 62 (2015) 21–29

magnitude of the variation due to the Swinson ﬂow was found to
be n = (0.1648 ± 0.0020)%
The radial density gradient of the galactic cosmic rays near the
Earth, based on our measurements may be estimated as described
below. The density gradient Gr ðPÞ of cosmic rays at a rigidity P may
be formulated as follows,jGr ðPÞj n ðPÞ
q sin v

ð23ÞHere, q is the gyro radius (Larmor radius) of the cosmic rays in the
IMF, and v the angle between the IMF and the direction of motion of
the Earth and is expected to be 45°. n ðPÞ represents the anisotropy
along the north–south direction perpendicular to the ecliptic plane
due to the Swinson ﬂow. As explained earlier, the median rigidity of
the cosmic rays responsible for the production of muons observedin the GRAPES-3 muon detector was P = 77 GV. From these parameters q was found to be 0.36 AU using the observed mean value ofthe IMF, B = 4.7 nT, measured near the Earth.
The radial density gradient of the cosmic rays at a rigidity of
77 GV was estimated to be 0.65% AU1, from these six years ofobservations with the GRAPES-3 experiment. The estimated density gradient from the GRAPES-3 data is comparable to the values
reported by other workers [2,9]. For example, the parametrization
given in [9] yielded a radial density gradient of 0.7% AU1 for arigidity of 77 GV. Similarly, the work in [2] implies a radial density
gradient of <1.0% AU1 at 100 GV.
5. Conclusion
The outward ﬂow of solar wind results in a gradient in the radial
density of galactic cosmic rays in the heliosphere due to the phenomenon of diffusion and convection. This gradient coupled with
the interplanetary magnetic ﬁeld results in a ﬂow of charged particles termed the Swinson ﬂow perpendicular to the ecliptic plane.
The Swinson ﬂow is manifested as a small anisotropy in siderealtime in the ﬂux of cosmic rays and its amplitude and phase were
measured using the GRAPES-3 tracking muon telescope data. The
large area (560 m2) of GRAPES-3 muon telescope resulted in thedetection of 7.5  1012 muons over the six year period during
2000–05. This high muon statistics was exploited to make a very
high precision measurement of the amplitude of the Swinson ﬂow,n = (0.0644 ± 0.0008)% which is a nearly 80r effect. The phase of
the maximum ﬂow was measured to be at a sidereal time of(17.70 ± 0.05) h. The phase of 17.7 h is close but different by
18 min from the expected value of 18 h. This difference of 18 minis small, yet is signiﬁcant at the level of 6r which is indicative of
the precision of the GRAPES-3 measurement. The mean interplanetary magnetic ﬁeld (IMF) was calculated to be 4.7 nT during thesix year interval 2000–05. The median rigidity of the primary protons responsible for producing the muons detected by GRAPES-3
experiment was 77 GV. From these values of Swinson ﬂow, mean

29

IMF and median proton rigidity the cosmic ray radial density gradient was estimated to be 0.65% AU1. This value is consistent withthe measurements reported by other groups.
Acknowledgment
We thank D.B. Arjunan, late S. Karthikeyan, S. Kingston, K.
Manjunath, S. Murugapandian, B. Rajesh, C. Ravindran, V. SanthoshKumar, R. Sureshkumar for their help in the operation, maintenance of the proportional counters and the associated electronics
for the GRAPES-3 tracking muon telescope. We thank G.P. Francis,
V. Jeyakumar, and K. Ramadass for their help in the repair andmaintenance of various mechanical components in the muon telescope. The easy availability of the OMNI data from the GSFC/SPDF
OMNIWeb interface at http://omniweb.gsfc.nasa.gov made it possible to determine the IMF direction which was a critical requirement of the data analysis. Use of the Kiel neutron monitor data
was important for our analysis for which we are grateful to the Kielgroup. This work was partially supported by the grants of the Solar
Terrestrial Environment Laboratory of Nagoya University, theInstitute for Cosmic Ray Research of University of Tokyo, Chubu
University, and the National Astronomical Observatory of Japan.
References[1] E.N. Parker, Phys. Rev. 110 (1958) 1445.
[2] D.L. Hall, M.L. Duldig, J.E. Humble, Space Sci. Rev. 78 (1996) 401.
[3] P.A. Isenberg, J.R. Jokipii, Astrophys. J. 234 (1979) 746.
[4] D.B. Swinson, J. Geophys. Res. 74 (1969) 5591;
D.B. Swinson, J. Geophys. Res. 81 (1976) 2075.
[5] K. Munakata et al., J. Geophys. Res. 104 (1999) 2511.
[6] F.J.M. Farley, J.R. Storey, Proc. Phys. Soc. A 67 (1954) 996.
[7] K. Nagashima, H. Ueno, Rep. Ionos. Space Res. Jpn. 25 (1971) 212.
[8] K. Nagashima, R. Tatsuoka, S. Matsuzaki, Nuovo Cimento C 6 (1983) 550;
R. Tatsuoka, K. Nagashima, Nuovo Cimento C 8 (1985) 320.
[9] S. Yasue, J. Geomagn. Geoelec. 32 (1980) 617.
[10] H. Moraal, Space Sci. Rev. 19 (1976) 845.
[11] W.I. Axford, Ap. J. Suppl. 90 (1994) 937.
[12] P.K. Mohanty et al., Astropart. Phys. 31 (2009) 24.
[13] P.K. Mohanty, S.R. Dugad, S.K. Gupta, Rev. Sci. Instrum. 83 (2012) 043301,
http://dx.doi.org/10.1063/1.3698089.
[14] S.K. Gupta et al., Nucl. Instrum. Methods A 540 (2005) 311.
[15] Y. Hayashi et al., Nucl. Instrum. Methods A 545 (2005) 643.
[16] S.K. Gupta et al., Nucl. Phys. B Proc. Suppl. 196 (2009) 153.
[17] H. Tanaka et al., J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201.
[18] T. Nonaka et al., Phys. Rev. D 74 (2006) 52003.
[19] P. Subramanian et al., Astron. Astrophys. 494 (2009) 1107.
[20] K.P. Arunbabu et al., Astron. Astrophys. 555 (2013) A139, http://dx.doi.org/
10.1051/0004-6361/201220830.
[21] C.C. Finlay et al., Geophys. J. Int. 183 (2010) 1216.
[22] <http://www-ik.fzk.de/corsika>.
[23] <http://134.245.132.179/kiel/main.htm>.
[24] W. Kolhörster, Phys. Z. 42 (1941) 55.
[25] H. Alfvén, K.G. Malmfors, Arkiv. Mat. Astron. Fysik. 29A (1943) 24.
[26] H. Elliot, D.W.N. Dolbear, J. Atmos. Terr. Phys. 1 (1951) 205.
[27] J.H. King, N.E. Papitashvili, J. Geophys. Res. 110 (2004) A02209, http://dx.doi.org/10.1029/2004JA010804.
<http://omniweb.gsfc.nasa.gov/html/
omni2_doc.html> .
[28] P.K. Mohanty et al., Pramana J. Phys. 81 (2013) 343.