Analysis of 25 mutual eclipses and occultations between the Galilean satellites observed from Brazil in 2009

The light curves of mutual eclipses and occultations between the natural satellites of a planet allow us to obtain high-precision position and relative motion from differential photometry, enough to detect weak orbital disturbing forces, such as tidal forces. The observations are made during the equinoxes of the planet. 
We studied 25 light curves observed in Brazil during the 2009 campaign of the Galilean satellites’ mutual phenomena. A narrow-band filter centred at 890 nm was used, strongly attenuating the Jupiter's scattered light. We fitted the occultation and eclipse light curves using semi-analytical models that take into account the gradual decrease of light over the shadow, the solar limb darkening and the solar phase angle. The Oren–Nayar reflexive model was used to map the inhomogeneous light scattering on the surface of the satellites. For the first time it is used in a work about mutual events. Here, we include the study that made us decide for this model. 
We measured the impact parameter, relative velocity and central instant with average precisions of 7.46 km (2.2 mas), 0.08 km s−1 (0.02 mas s−1) and 0.42 s (6.13 km), respectively. The fit precision of the normalized light-curve fluxes ranged between 0.4 and 4.4 per cent.

In this context, the photometric observations of mutual phenomena have an important advantage over traditional astrometric observations. Although such observations can only be done from time to time, during the equinox of the planet, they provide much better precisions. Indeed, for the Galilean satellites, the precision obtained with the current astrometric techniques, for a single satellite position, ranges between 134 and 170 mas (milliarcsecond) (Kiseleva et al. 2008). For relative distances between satellite pairs, the 30 mas error level can be achieved (Peng et al. 2012). On the other hand, for past mutual phenomena, the errors can easily reach the 20-30 mas range, and many times not exceed 5 mas (Emelyanov 2009). In this work, the average precision was even better, 3.10 mas.
Mutual phenomena between the Jovian satellites were observed for the first time in 1973 (Aksnes & Franklin 1976), and again in 1979, followed by the observation of the mutual events of the Saturnian satellites in 1979-1980(Aksnes et al. 1984. A detailed history of all the past campaigns of mutual events between the satellites of the giant planets, as well as a compilation of past reduction methods and recent improvements in the fit of light curves can be seen in Emelyanov (2009) and in the references therein.
Because of the importance of mutual phenomena, many researchers have worked in their prediction and organized observational campaigns to study both of these systems (see Aksnes & Franklin 1990;Thuillot & Arlot 1996;Arlot, Lainey & Thuillot 2006; for the more recent works). In 2007, mutual phenomena between the Uranian satellites were observed for the first time with CCD technology (since they occur every 42 yr), and provided important astrometric results (see, for example Birlan et al. 2008;Assafin et al. 2009 andChristou et al. 2009).
Based on the prediction of the events between the Galilean satellites for 2009 (Arlot 2008), a Brazilian campaign was organized, involving researchers from four institutes. Observations were carried out at the Observatório do Pico dos Dias (OPD), managed by the Laboratório Nacional de Astrofísica (LNA), Itajubá, Brazil (IAU code 874). We successfully observed and analysed 25 events (13 occultations and 12 eclipses). Here, we present the observational procedures, the data reduction and the analysis of these events.
In Section 2, we describe the campaign's programme and observations. In Section 3, we present the photometry and the lightcurve fitting procedure, with a complete description of the analytical models used. We give the results in Section 4. Section 5 brings a discussion about the available light reflectance models, and the introduction of the Oren-Nayar's model, adopted in this work. Conclusions are set in Section 6. In Appendices A and B, we detail the calculations involved in the light-curve fits, related to the adopted analytical models. In Appendix C, we display the fits to the entire set of 25 light curves analysed in this work.

Programme
The observational campaign in Brazil was the result of a collaboration between four national institutes. The predicted events for 2009-2010 (Arlot 2008) were available at the portal of the 'Institut de Mécanique Céleste et de Calcul desÉphémérides' (IMCCE) 1 from the Observatoire de Paris. We selected all the visible events for the OPD/LNA observatory 2 with a predicted non-zero flux drop and elevation above 30 degrees. This resulted in 50 events distributed in 45 nights spread over nine months. We attempted observations for all these events.

Observations
We lost 23 events due to bad weather conditions. Also, from the remaining 27 mutual phenomena successfully observed in 23 nights, two were quasi-simultaneous, superposed in the same light curve. These two events will not be studied in this work. The resulting 25 events analysed in this work are divided into 12 eclipses and 13 occultations. The satellites Io, Europa and Ganymede were all eventually observed in an occultation or in an eclipse and Callisto was observed only in eclipses. 1 ftp://ftp.imcce.fr/pub/ephem/satel/phemu09/phemu09_132ts.txt 2 ftp://ftp.imcce.fr/pub/ephem/satel/phemu09/visibility/vtri-itajuba.txt The events were observed at the OPD/LNA (λ = +45 0 32 57 , φ = −22 0 32 22 , h =1860 m, IAU code = 874). The observations were carried out using the 0.6 m diameter Zeiss telescope, f/12.5. For one night (20/06), the 1.6 m diameter Perkim-Elmer telescope, f/10, was used.
Owing to methane in the Jupiter's atmosphere, the planet has a strong absorption of light between 880 and 900 nm, which causes the planet albedo to drop to 0.1 in this region (Karkoschka 1994(Karkoschka , 1998. As it does not occur with the satellites, a narrow-band filter at these wavelengths was used to strongly minimize the scattered light of Jupiter. The effect is shown in Fig. 1. The planet and its satellites present about the same brightness. This 'methane' filter is centred at 890 nm with a bandwidth of 20 nm. We used two back illuminated CCD detectors. The EEVCCD detector (model 02-06-1-206) with 385 × 578 square pixels of 22 μm, hereafter CCD 301, was used in 23 out of the 25 events studied in this work. The EDVCCD detector (model 47-20) with 1024 × 1024 pixels of 13.5 μm was used in the two events observed at the Perkim-Elmer telescope (2006CeI-1 and 2006CeI-2, see Table 1). Since mutual phenomena are short-term events (typically a few minutes long) with relative satellite velocities above 6 km s −1 , they demand short exposures for achieving a time resolution corresponding to a spatial resolution of a few kilometres. In practice, depending on the specific relative speed and on the weather conditions, the exposure times ranged between 1 and 3 s. This granted satellite's ADU (analogic-to-digital unit) peak counts of about half the CCD maximum, allowing for high signal-to-noise ratios (S/N) at the optimal linear count ranges of the CCDs. The electronics of the detectors were set to do simultaneous integration and charge transfer (frame-transfer mode), eliminating the readout overhead between acquisitions.
Observations were made in order to always keep the same photometric calibrator in the field of view (FOV). In eclipses, sometimes the eclipsing satellite was the calibrator. In the absence of suitable satellites, Jupiter -and sometimes a spot on its surface -was used instead. Yet, due to the large separation of the targets, two events were observed without a calibrator: an occultation of Io by Europa on August 07 and an occultation of Ganymede by Europa on August 12.  ), and the event type ('e' and 'o' stand for eclipse and occultation, respectively). Also indicated are the objects used as calibrators (Cal.) in the differential photometry: J stands for Jupiter, S for a spot on Jupiter and N means no calibrator available. We also give the seeing, the total number of images used, the zenith distance z and the solar phase angle. There was no prediction in (Arlot 2008) for the event 1208GoE.
The FOV was 4 arcmin × 4 arcmin, with a pixel of 0.6 arcsec size. Seeing was typically in the range 1-2 arcsec.
Observations started 30 min before and ended 30 min after the predicted instants provided, respectively, for the start and end of the event. This procedure aimed to obtain well resolved images of the satellites involved, with enough angular separation to measure their individual fluxes, as close as possible to the event, to determine the ratio of albedos (see discussion in Section 3.2). Table 1 gives information for the 25 events analysed in this work, indicating the targets, calibration object, seeing, zenith distance, number of images and solar phase angle.

P H OTO M E T RY, L I G H T-C U RV E F I T T I N G
In this section, we will describe: (a) the differential aperture photometry applied to generate the observed light curves; (b) the technique used for determining the ratio of albedos; (c) the analytic models and numerical computations utilized to fit the observed light curves.

Photometry
All images were corrected by bias and flat-field using the IRAF package (Butcher & Stevens 1981). Differential aperture photometry was performed using the PRAIA package described in Assafin et al. (2009). The flux of the target was measured relative to a calibrator, usually another satellite or satellites, sometimes the eclipsing one (in eclipses). The resulting measured target/calibrator flux ratio is practically free from the sky transparency variations owing to anomalous atmospheric extinction, as this affects all objects in almost the same manner in the FOV. The S/N, the flux ratio error and the seeing are calculated and stored together with the mid-time instant of the measurements. After measuring the flux ratio for all the images, we obtained the observed light curve of the event.
In order to maximize the S/N, for each event we tested the optimum size of the aperture radius for the flux determination of each object, and for the sky background annulus. We varied their sizes for each series of observations, and plotted the resulting sets of provisional light curves. We choose the aperture values and sky background annulus which resulted in the light curve with the least flux dispersion. Typically, the sky background annulus around the objects had an internal radius 5 pixels larger than that of the aperture circle, and a width of 5 pixels. We did not verify a clear linear relation between the best aperture radius size and the seeing of the night for the events of this campaign.
From Table 1, we notice that for some events, Jupiter or a spot on its surface were used as calibrator. We verified the validity of such calibrators by comparing the light curve of an event that had a satellite as calibrator, with the light curve of the same event obtained with Jupiter or a spot as the calibration object. Although we add a certain amount of noise in the flux ratio, the procedure proved to be effective.
In two cases, no calibrator at all was available in the FOV. For these events, we fitted a polynomial to the light curve outside the flux drop, and flatted the entire curve to eliminate the continuous variation of the flux, due to the atmospheric extinction. The procedure proved to be satisfactory, since the sky conditions were very stable. In both cases, the light curves were thus directly obtained from the fluxes of the targets.
As a final step, all observed light curves were normalized by fitting a polynomial curve of nth degree (usually n = 1) outside the flux drop, so that the flux ratio was set to 1 outside the events.

Ratio of albedos
The ratio of albedos is not relevant for eclipses when the satellites' photometric fluxes are separately measured, as is our case. Therefore, the following discussion concerns only to occultations.
The ratio of albedos is a physical parameter that is related to the reflectivity of the satellites. Since the albedo is the ratio between the received and reflected light by the satellites, the ratio of albedos has great influence on the flux ratio shape of the light curve of an occultation. In particular, it has a high correlation with the impact parameter, which is the minimum distance at apparent closest approach. Because of this correlation, it is highly recommended that we determine the ratio of albedos independently from the reduction process, instead of trying to fit it from the light curve, together with the other parameters.
There are maps of albedos for the Galilean satellites provided by the Voyager and Galileo probes (see, for example, Vasundhara et al. 2003). But they were not generated in the same bandpass (filter) as our observations. And there is no reliable process to convert these albedos to our methane bandpass, or to other wavelengths, because of photoclinometry effects, such as shadows cast by the relief at the time of the probe's measurements, which are not present in the ground-based observations, and vice-versa.
Instead, using the procedure first mentioned in Assafin et al. (2009) (see Section 5), also followed by Emelyanov et al. (2011), we determined the ratio of albedos directly from the observations, using images before and after the occultations, with the involved satellites fully resolved in the images. Since the albedo is the ratio between the satellite's photometric flux and area, we can write the ratio of albedos of satellites 1 (A 1 ) and 2 (A 2 ) as: where F 1 and F 2 are the photometric flux of occulted and occulting satellites, respectively, obtained from the images taken near the event, and S 1 and S 2 are their areas. We did not use the generalized form of equation (1) proposed by Emelyanov et al. (2011), because it increases the error in the determination of the albedo ratios, by (unnecessarily) introducing one more flux quantity (the sum of satellite fluxes F 12 or F 21 measured together), which in turn cannot be measured without avoiding large contributions of sky background and Poisson noise.

Light-curve fitting procedure
We rigorously reproduced the geometry of the mutual events by formulating theoretical models and performing numerical computations. We then could determine the apparent configuration of the satellites and of the Sun (including shadows), and thus be able to infer the amount of light received by the observer at any instant. This flux in turn was compared and fitted to the observed light curves in an iterative procedure. As a result, we determined the observed apparent relative orbital paths of the target satellites, described in terms of known parameters: the impact parameter (the minimum apparent distance between the geometric centre of the satellites, at the apparent closest approach in the sky plane), the relative velocity (tangent to the point associated with the impact parameter, measured in the sky plane) and the central instant (instant of time associated with the impact parameter). It is important to emphasize that the impact parameter is defined by using the geometric centre of the satellites which, unlike the photocentre, do not depend on the solar phase angle. The satellites' ephemeris can also be represented by these parameters, so that we could compare the results from the observations with the current established ephemeris for the Galilean satellites.
We assume an apparent linear path for the relative motion between the target satellites, during the short time interval of the events. We model the light curves of the events based on a rigorous formulation of the geometry of the relative orbits and of the effects of light reflectance over the satellites' surfaces. The normalized flux ratio cannot be directly expressed in terms of a simple analytical function. For each instant of observation, we calculate the flux by numerically simulating the two-dimensional luminosity profile of the satellites' discs, projected on to the sky plane perpendicular to the line of sight. These profiles depend on the type of the event (occultation or eclipse), and on the positions of the satellites and of the Sun. As a result, we derive a simulated light curve, which works as the model. This light curve is parametrized by the relative orbital parameters described above. The observed light curve is then fitted to this simulated light-curve model by an iterative non-linear least-squares procedure, which follows the Levenberg-Marquardt method. The derivatives of the simulated light-curve model with respect to the orbital parameters are numerically computed by varying the parameters. In the process, as the flux of the simulated lightcurve model converges to the flux of the observed light curve, so do the values of the orbital parameters being adjusted. After a convergence of 1 per cent is reached in the chi-square (flux) between the simulated and observed light curves, we obtain the desired orbital parameters. The spatial resolution of the luminosity profile is set (limited) by the time resolution and photometric measuring errors of the observations, thus at the same time avoiding loss of precision, and optimizing computation speed. Typically, we used resolutions corresponding to 10 × 10 km areas in the sky plane.
We fit the projected relative velocity, the central instant, and the impact parameter with this procedure. Initial guess values for these parameters are necessary to start the procedure. They are solely obtained from the observed light curve itself, which makes the process ephemeris independent. On the other hand, we use the formalism above to rigorously compute relative velocity, central instant and impact parameter, based on a given ephemeris. This translation allows for a direct and unbiased comparison between the ephemeris-based and the fitted parameters, which model the observed light curves of the mutual events (see Section 3.4).
The initial guess values are calculated from the observed light curve as follows. The central instant is obtained from the instant associated with the observation with minimum flux. The relative velocity is computed by considering the diameter of the occulted/eclipsed body and the time interval from ingress to egress, estimated from the observed light curve. The impact parameter is derived in an iterative process. We increment the impact parameter values from zero to a given limit, which is the sum of both satellites' radius for occultations, or the penumbra size plus the radius of the eclipsing satellite, in the case of eclipses. For each incremented value, we compute the simulated flux associated with the estimated central instant and to four other points, uniformly located between ingress and egress. Then, we store the chi-square, which is the square root of the ratio between the sum of the differences between the computed and observed fluxes, and the degrees of freedom (N points minus the number of parameters). We keep the impact parameter value corresponding to the least chi-square found.
In the following, we detail the modelling and calculations necessary for the numerical simulation of the light curves, separately carried out for occultations and eclipses.

Light-curve simulation for occultations
Three factors have great importance on the theoretical modelling of occultations: the ratio of albedos, the reflectance law and the solar phase angle. The first two determine the luminosity intensity profile of each point in the projection, while the third one shapes the apparent discs at the sky plane.
We adopt the simplified version of the model proposed by Oren & Nayar (1994), which consists of a generalization of Lambert's law, with the surface roughness seen as a set of facets with different inclinations. We use the simplified version because, according to the authors, it describes well the problem and, for the purposes of this work and achieved precisions, there was no need to consume the additional processing time, required to use the complete version.
The Oren-Nayar's model has been used in works on threedimensional reflective surfaces since its formulation (see Yinlong 2007). This model is close to the actual reflection profile of an illuminated curved surface, providing a satisfactory fit to the light curves. To our knowledge, it is the first time that this model is used in a work about occultations and eclipses in mutual events.
We must know the geometric orientations of the incident and reflected (emergent) rays, respectively L i and L e , which are determined by four angles, as shown in Fig. 2. Here, ψ i and ψ e represent the angles between the normal vector η 1 to the surface, respectively, Downloaded from https://academic.oup.com/mnras/article-abstract/432/1/225/1117425 by guest on 25 July 2018 Figure 2. The spacial orientation of the incident and reflected (emergent) radiances, geometrically represented by the light rays L i and L e , for a point in the satellite's surface. ψ i and ψ e represent the angles between the normal vector to the surface, η 1 , respectively, with L i and L e . φ i and φ e are obtained by projecting the rays at the plane tangent to the satellite's surface and taking the angles between these projections and the vector η 2 , parallel to the tangent plane and in the north-south direction (see Appendix A and Fig. A1).
with L i and L e . φ i and φ e are obtained by projecting the radiances at the plane tangent to the satellite's surface, and taking the angles between these projections and the vector (η 2 ), parallel to the tangent plane and in the north-south direction (see Appendix A and Fig. A1).
Thus, considering the albedo A of the body, and given the incoming and reflected radiances, geometrically represented by the light rays L i and L e , according to this model we have where The surface roughness is represented by the variable σ , which is the aperture angle between the facets. It ranges between 0 (perfectly smooth surface) and π/2 (entirely rough surface). Here, we used σ = π/2. In Section 5, we discuss the choice of the reflectance model, and the choice of σ in the fitting of the light curves.
The Oren-Nayar's model requires the determination of the spatial orientation of the incident and reflected radiances stated in equation (2). This requires the computation of the angles shown in Fig. 2.
To determine these angles, we used state vectors V , consisting of the topocentric and heliocentric position and velocity of the satellites and of the Sun (see Fig. 3). They were generated using the DE418, combined with the file NOE-5-2010-GAL.a.bsp through the SPICE information system (Acton 1996). This file, available on is the satellite 2 heliocentric position vector. For clarity, the velocity components were omitted. the web, 3 provided the IMCCE's ephemeris (theory NOE-5-2010) used in this work for the four Galilean satellites (Lainey et al. 2009). Note that the errors in the state vectors from the current ephemeris have a negligible propagation effect in the computed quantities, and can be virtually ignored.
From Figs 2 and 3, we have L ij = V S j and L ej = V j , where j is 1 for the occulted satellite (satellite 1), and 2 for the occulting satellite (satellite 2). Thus, from the cross product and dot product between the state vectors, we have We show in Appendix A how to compute the normal and tangent vectors η 1 and η 2 . The topocentric and heliocentric state vectors provide, for each instant, the direction of the incident and scattered light rays (radiances). The orientation of the vectors η 1 and η 2 depends of the position of the point in the satellite's surface.
We use these expressions with the Oren-Nayar model to determine the luminous intensity of each point of the disc projected in the sky plane. Note that the effect of the solar phase angle in the projection is naturally accounted for by the reflectance model. Also, the procedure easily identifies when a point in the disc of the occulted satellite lays behind a point in the occulting satellite in front of it. The integration of the (normalized) flux over the two discs gives the total simulated flux of the target satellites for a given instant. We then have the flux of a single point in the simulated light curve.
We finally obtain the entire simulated light curve by repeating the procedure for all the instants corresponding to an actual observation.
Importantly, all geometry involving state vectors takes place in a timeless space. Thus, the instants used to generate the state vectors of each satellite and the sun are the instants of each image corrected by the light travel time for each position.

Light-curve simulation for eclipses
The shadow of an eclipse is composed by two regions, the umbra and the penumbra. The flux is zero for points inside the umbra. Inside the penumbra, the flux is attenuated, as compared to a 'no-eclipse' situation. Therefore, it is fundamental to rigorously determine the umbra and penumbra regions, for a very precise fitting to the observed light curve of an eclipse.
We geometrically determined the radius of the umbra and penumbra, and semi-analytically derived the gradual decrease of light along the penumbra.
The radii of the umbra (R U ) and penumbra (R P ) at the path of the eclipsed satellite were determined from the state vectors, using basic geometry (see Fig. 4): where R S and R 2 are the radii of the Sun and of the eclipsing satellite, respectively. = |V S1 |, so that the angle between the heliocentric vectors V S1 and V S2 (see Fig. 3) is arccos V S1 ·V S2 |V S1 ·V S2 | . Also, O S E = V S1 ·V S2 |V S2 | . The solar plane along AO S B contains the Sun, and is normal to the heliocentric vector of the eclipsing satellite. See also the discussion in Appendix B. We obtained the flux profile of the penumbra by determining the fraction of extinct sun's light along the path of the eclipsed satellite. This was made in two steps.
First, we consider a fictitious observer placed at a point in the penumbra looking towards the Sun, and analytically determine the fraction of the disc of the Sun, which is covered by the disc of the eclipsing satellite. For that, we used the expressions from Assafin et al. (2009), which determine the common area (S c ) between two discs of radii R S and R SS overlapping each other, as a function of the distance d between their centres (see Fig. 5): where Here, R S (the Sun radius) and R SS are taken at the solar plane, which is the plane containing the Sun, perpendicular to the heliocentric vector of the eclipsing satellite (see Fig. 4). R SS is the radius of the eclipsing satellite, projected at this plane from the fictitious observer point of view at the penumbra. The calculation of R SS , and the use of R SS , R U and R P in the determination of the actual eclipse shadow cast in the observation plane, including the effect of solar phase, are explained in Appendix B.
The second step consists in taking into account the solar limb darkening, at the uncovered fraction of the solar disc, to calculate the amount of light that is actually received by the satellite. For that, the uncovered solar disc radial profile is numerically represented by rings. Each ring has a light intensity value coherent with the darkening profile proposed by Hestroffer & Magnan (1998). The number of rings was chosen so that the light intensity difference between two consecutive rings was less than 1 per cent, resulting in a smooth limb darkening profile.
Thus, from the obtained penumbra light profile, we determined the incoming and outcoming fluxes L i and L e with the reflectance model, for any point of the eclipsed satellite in the sky plane.
The points of the eclipsed satellite disc projected in the sky, bathed by the umbra, have zero flux. The projected satellite points inside the penumbra have intermediate fluxes, computed by the procedure explained in this section and in Appendix B. Those satellite points outside the penumbra (if any) have fluxes computed by the same procedure described in the previous section. In this way, we derive the luminosity intensity at each point of the satellite disc projected in the sky plane. Integration of all fluxes furnish the (normalized) simulated flux for a given instant. Similarly as before, after considering all the instants, we obtain the entire simulated light curve for the eclipse.

Ephemeris offsets and ephemeris-based light curves
From state vectors, applying basic plane geometry, we can obtain the impact parameter, the central instant and the relative velocity between two given satellites. Since the state vectors are generated from the ephemeris, these values may be regarded as ephemeris based, and can be compared to the values obtained from actual fits to the observed light curves. A simple transformation translates the differences found into right ascension and declination 'observed minus ephemeris' offsets. In fact, this information is stored in the fitting procedures.
The fitting software constructs a light curve from values for the impact parameter, relative velocity and central instant. This allows for a visual inspection during the fitting procedure, from starting values to the final results themselves. The ephemeris-based computed values are also used to generate ephemeris-based light curves, that can be compared to the observed and fitted ones.
All this allows for both a qualitative and quantitative comparison between observation and ephemeris.

R E S U LT S
We fitted the light curves of 25 events and obtained the values and error estimates for the impact parameter, the relative velocity and the central instant, following the procedures described in Section 3. We also calculated these parameters from the ephemeris, as explained in Section 3.4. This allowed us to compare the expected results from the ephemeris with the observed ones. Figs 6 and 7 illustrate two examples of adjustment of light curves from an eclipse and from an occultation. The ephemeris-based light curves are also displayed. In Appendix C, we illustrate all the events. Figs C1 to C3 (occultations) and C4 to C6 (eclipses) sample the events according to the quality of the observations and adjustments, i.e. 'excellent', 'good' and 'regular'.
We show in Table 2 (for occultations) and in Table 3 (for eclipses) the parameters and respective errors obtained from the light-curve fits of all the events. For the impact parameter and the velocity, values are listed in kilometres and kilometres per second, and in milliarcseconds and milliarcseconds per second, respectively. For the central instant, in UTC (Universal Time Coordinated), the label of each event indicates the day and month. We list the instant in hours, minutes and seconds, and the error in seconds and in kilometres (by the use of the relative velocity in kilometres per second). We also list in Tables 2 and 3 the photometric error, based on the dispersion of light-curve points outside the events, and the mean error in flux ratio, computed from the light-curve fits. The ratio of albedos (and errors) used for the light-curve reduction of the occultations (see Section 3.2) are also given in Table 2.
We compared the results with the ephemeris published by the IMCCE, currently considered the most accurate representative for the Jovian system. We list in Tables 2 and 3 the differences between the fitted and the ephemeris-based parameter values. We also list the ( αcos δ, δ) orbital offsets in the sense 'observation − ephemeris'.

T H E R E F L E C TA N C E L AW
The choice of the reflectance law has great influence in the lightcurve fit of the mutual events. This is because the satellite's photocentre is significantly offset relative to its geometric centre, if the solar phase angle is not negligible. Therefore, the choice of the reflectance model affects the relationship between the impact parameter and the light-curve's depth. Another possible side effect is the bad determination of the central instant. Thus, the choice of a simple or incorrect reflectance model may lead to a less precise, possibly inaccurate determination of the parameters, depending on the photometric quality of the light curve. As a consequence, the mean error of the fits in flux ratio may increase. We used the computed mean errors of some selected light curves as a guide to measure the adequateness of some tested reflectance models. The final choice was made according to the quality of the fit obtained from these light curves.
Knowing that the Poisson noise is proportional to the square root of the intensity of the object's light flux, a decrease in brightness, as observed close to the central instant of some events (mostly central eclipses), causes a decrease in the noise, and hence, in the dispersion of the points near the bottom of the light curve. A narrowing is observed in the light curve near the central instant for such events. Fig. 8 displays an example. Note. The results are arranged in two rows for each event. The values in parentheses are their respective errors. The impact parameter and the relative velocity are listed in the first line, respectively, in kilometres and kilometres per second, and in the second line in milliarcseconds and milliarcseconds per second. For the central instant, in UTC, the label of each event indicates the day and month, and the first line gives the instant in hours, minutes and seconds. The second line (σ MI ) lists the mid-time instant error in kilometres, by the use of the relative velocity in kilometres per second. IP, RV and MI are the respective differences between the fitted parameters and the ephemeris-based values in the sense 'observation − ephemeris'. In columns 8 and 9, we also list the ( αcos δ, δ) orbital offsets in the sense 'observation − ephemeris'. The photometric errors are based in the dispersion of the light-curve points outside the events. The σ fit values are the mean error in the flux ratio, computed from the light-curve fits. The listed ratio of albedos were computed following the procedures described in Section 3.2. For 1208GoE, the ratio of albedos was not directly determined from observations before and after the event. Instead, we computed it by multiplying the values of the ratio of albedos of the event 2704GoI, and the average value for Io/Europa ratio of albedos obtained from the events 0708IoE, 2208IoE and 1609IoE. These events were observed close to the 1208GoE occultation.
In the Poisson noise regime (as is the case of our observations), we should expect that the photometric error is minimum at the bottom of the light curve, thus making this part of the light curves suitable for the analysis of the reflectance model. If the reflectance model does not provide a good fit, the bottom part of the theoretical light curve should be slightly offset from the observed one, due to the bad determination of the impact parameter and, possibly, also of the central instant. This effect is shown in Fig. 9. Therefore, besides the inspection of the mean error of the flux ratio, computed from the light-curve fits, this behaviour is another factor that helps in the choice of the reflectance model.
The most basic and commonly used light scattering model is the simple Lambert's scattering law. It assumes an infinite-sized source of light, and a reflecting surface that scatters the incident light equally in all directions. Therefore, the satellite's surface luminosity can be represented by a homogeneous grey disc. In this case, one way to account for the solar phase angle is to describe the format of the apparent disc as a combination of a semi-circle with a semi-ellipse. However, the high photometric quality of the observed events makes this model too simplistic and, therefore, inappropriate to provide a high-quality fit to our light curves.
In this work, we considered the most complex and complete reflectance models used in the literature, that take into account the direction of the incident light for a finite-sized source, for example. First, we studied the generalized version of Lommel-Seelinger's law, where the ratio between the radiance of the incident light (L i ) and the observed object's radiance (L e ) is given by (Buratti & Veverka 1984) Here, f(α) = 1 + Bα + Cα 2 is the phase function of the surface, where α is the phase angle, ψ i and ψ e are the same as in Fig. 2, A is a function of α and the parameters B and C depend of the satellites surface features. Unfortunately, in the literature, there was Note. Read definitions in Table 2. We do not list the ratio of albedos, as this is not used in the reductions of eclipses. no information about these parameters for the Galilean satellites (except for Europa), for the wavelength range covered by this study. Therefore, this model could not be further tested. We then tested the non-generalized version of Lommel-Seelinger's law (equation 15) (Aksnes, Franklin & Magnusson 1986), in an attempt to eliminate the need of information about the features of the satellites. The law can be described as cos ψ i cos ψ e cos ψ i + cos ψ e ds (15) where, d(L e /L i ) is the amount of radiance scattered from a surface element ds, ψ i and ψ e are the same as in Figs 2 and B is a normalization constant that, for this work, is assumed to be the value of the ratio of albedos, aiming to maintain the aspect ratio between the satellites fluxes. This alternative did not satisfactorily solve the problem of the poor determination of the central instant, evidenced by Fig. 9 and, therefore, a different model was needed.
Hapke's scattering law (Hapke 1981(Hapke , 1984) is a sophisticated model that has been used in some works. Unfortunately, we could not apply it, due to the severe lack of information on many of the parameters of that law, for the wavelength band of our observations -as was the case with the generalized version of Lommel-Seelinger's law. We notice, however, that Emelyanov (2009) (see Section 4, table 2 in that article), reports that both the Hapke's and Lommel-Seelinger's laws have a similar performance.
We finally found the solution as a generalization of Lambert's law (Oren & Nayar 1994), that takes into account not only the direction of the radiances, like the other models, but also the surface roughness. The roughness is represented by a factor of σ , which ranges between 0 and π/2, with π/2 equivalent to the surface roughness of a full Moon. The Oren-Nayar's model proved to be efficient in describing the profile of spherical surfaces illuminated by a finite source, and then we decided to test it.
The first test was to verify the model's capacity to solve the problem pointed out in Fig. 9. The results were highly satisfactory, as shown in Fig. 10.
We then verified if there were a strong dependence of the model with the parameter σ describing the surface roughness. If we found a high dependence, we should consider fitting this parameter too. However, changing the value of σ did not produce significant changes in the parameters of the light curve. Also, no significant changes occurred in the mean error of the flux ratio (σ fit ). This indicates that the Oren-Nayar's model is a very robust reflectance model, even for high-precision photometric light curves. One advantage of this model is that no ad hoc albedo-wavelength-related parameters need to be set for the satellites, as is the case of other reflectance laws used in the literature. Thus, the Oren-Nayar's model is a very interesting alternative, that can be used in the light-curve fit of mutual events of any giant planet, without a priori photometric knowledge in the satellites. Taking into account the recommendations in Oren & Nayar (1994), we used σ = π/2 in our fits. The simplified version of the model was used. Due to our photometric precision, it furnishes the same results as the complete model, in much less computing time.

C O N C L U S I O N S
We presented the results obtained from the observation, reduction and analysis of 25 mutual events registered during the Brazilian campaign for the observation of mutual phenomena between the Galilean satellites.
The narrow-band filter at 890 nm, combined with the differential photometry and the refined procedures in the fitting of the light curves developed in this work, resulted in average precisions of 80.1 m s −1 (0.023 mas s −1 ) and 7.46 km (2.19 mas) for the relative velocity and impact parameter (relative positions), respectively, and 0.42 s (6.13 km) for the central instant.
The extensive use of numerical procedures with analytical and semi-analytical models, allowed for a complete, rigorous implementation of the complex geometry involved in describing the flux drop in occultations and eclipses. The modelling of the solar limb darkening and the implementation of the computation of the gradual decrease of light over the shadow in an eclipse are examples. This also allowed for the analysis of different reflective laws and models in a fast and direct way. Through this analysis, we were able to highlight and study the relation between the impact parameter and the reflective model. From this study, we could finally adopt a generalized reflectance model -the Oren-Nayar's model, used for the first time in a work on occultations and eclipses in mutual events. This model is well suited for fitting high-precision light curves, and does not depend on wavelength and other a priori photometric knowledge of the satellites, contrary to some reflectance models currently in use. We emphasize that our developed lightcurve fitting procedures are ephemeris-independent, thus allowing for an unbiased comparison with the current orbit implementations available for the Galilean satellites.
The tools and techniques developed and used in this work, allowed for a comprehensive analysis of the many factors inherent to the nature of these events. They will be of great utility in future campaigns, not restricted to the Jupiter system. The use of more powerful detectors will possibly make it easier to study new refinements in the analysis.

A P P E N D I X A : C O M P U T I N G V E C TO R S η 1 A N D η 2 F O R S I M U L AT I N G L I G H T C U RV E S
In the occultation case, and outside the umbra and penumbra in the eclipse case, we sweep the two-dimensional sky plane, integrating the reflected (emergent) radiance L e , in order to determine the flux value at each (x, z) point in that plane. In this way, we obtain the light intensity profile of the satellite discs projected in the sky plane -a necessary step in the simulation of the light curves. We obtain L e by equation (2), following the Oren-Nayar's reflectance model. The angles in equation (2) are computed from equations (6) to 9, which depend on the normal and tangent vectors η 1 and η 2 . Here, we show how we obtained η 1 and η 2 using spherical geometry.
Let us consider a spherical coordinate system with origin at the centre of the satellite (Fig. A1). The xz plane of the associated Cartesian coordinate system coincides with the sky plane (Fig. A2). Thus, a point (x, y, z) in the spherical surface of the satellite has spherical coordinates: where R is the satellite's radius. The Cartesian components of the vectors η 1 (r) and η 2 (θ ) are, thus: η 2 = θ = cos θ cos ϕx + cos θ sin ϕŷ + sin θẑ.
Now, we conveniently write η 1 and η 2 in topocentric coordinates, using the state vectors defined in Section 3.3.1. First, we define the x-axis parallel to the projection of the relative velocity vector on the sky plane. Then, we set the y-axis antiparallel to the satellite's  position vector, to ensure the perpendicularity to the sky plane (see Fig. A2). Finally, the z-axis results from the cross product between them. Thus, for each satellite j (j = 1, 2), we have where V Rel is the relative velocity vector, obtained from the difference between the topocentric velocity vectors of the satellites.
Equations (A6) and (A7) allow for computing the (x, y, z) Cartesian components of the vectors η 1 and η 2 , as a function of topocentric state vectors, and by inputting (x, z) values. This allows for sweeping the (x, z) coordinates of the satellite's circular discs projected in the sky plane. This sweeping is done by choosing a (x, z) resolution for the construction of the light profile of the discs, in a compromise between computational efficiency and accuracy.

A P P E N D I X B : C O M P U T I N G T H E F R AC T I O N O F S U N L I G H T I N E C L I P S E S
We compute the fraction of sunlight in eclipses by solving equation (12). For that, we must determine the radius R SS of the eclipsing satellite disc projected in the solar plane, for a fictitious observer placed at a point within the penumbra. We also need to calculate the distance d between the centres of this disc and the Sun's disc. For that, we use the shadow cones represented by the triangles AH p EBH p and IDJ (see Fig. 4; see also the discussion at the end of this appendix). Thus, we have where R CP is the distance of the fictitious observer to the umbra centre, along the path of the eclipsed satellite. It is convenient to express this distance as a function of the corresponding projected (x, y) coordinates in the observation plane, which is the sky plane perpendicular to the topocentric vector and containing the umbra centre. For that, it is necessary to consider two key angles. One is the solar phase angle between this plane and the solar plane. The other is , the angle between the observation plane coordinate system (x, y) and the relative orbital plane of satellites containing V Rel (see Fig. B1). Thus, in terms of these observation plane quantities, R CP can be expressed as R CP = [x 2 (cos 2 cos 2 + sin 2 ) + y 2 (sin 2 cos 2 + cos 2 )] 1/2 (B3) where, = |V 1 × V S1 | |V 1 ||V S1 | and = |(V S1 × V 1 )(V Rel × V 2 )| |V S1 × V 1 ||V Rel × V 2 | .
The (x, y) are the coordinates of a point in the penumbra with respect to the umbra centre at the observation plane. Using equations Figure B1. Coordinate system at the observation plane where the y-axis is normal to the plane containing the solar phase angle.
(B1) to B4, we sweep the (x, y) coordinates of this plane, considering the regions delimited by the computed umbra and penumbra radii, R U and R P . Taking the solar limb darkening into account, we then determine the effective light L i received by the eclipsed satellite, for each point in the sky. We then use the reflectance model to obtain the reflected light L e .
Comparing the shadow cones represented by the triangles AH p EBH p and IDJ in Fig. 4, we note that, for a fictitious observer located inside the umbra and penumbra (point D), the radius R SS of the disc projected in the solar plane undergoes a small variation, as we move along the segment EG. The amount of this variation depends on the distance between the satellites and the Sun, and on their radii. For all the events of this work, the maximum variation of R SS remained below 0.15 per cent, after considering the two extreme cases: observer located at the centre of the umbra (point E), and located at the edge of the penumbra (point G). This negligible effect was thus ignored, and we adopted the simple calculation of R SS using the triangle AH p EBH p alone (see Fig. 4).

A P P E N D I X C : L I G H T C U RV E S O F T H E E V E N T S
Here, we present the reduced light curves of the events treated in this work. The figures are divided by the quality of the light curves and respective events' codes are indicated in each figure.    This paper has been typeset from a T E X/L A T E X file prepared by the author.