# Accretion Disk Size Measurement and Time Delays in the Lensed Quasar WFI 2033–4723

###### Abstract

We present 13 seasons of -band photometry of the quadruply-lensed quasar WFI 2033-4723 from the 1.3m SMARTS telescope at CTIO and the 1.2m Euler Swiss Telescope at La Silla, in which we detect microlensing variability of mags on a timescale of 6 years. Using a Bayesian Monte Carlo technique, we analyze the microlensing signal to obtain a measurement of the size of this system’s accretion disk of at , assuming a inclination angle. We confirm previous measurements of the BC and AB time delays, and we obtain a tentative measurement of the delay between the closely spaced A1 and A2 images of days. We conclude with an update to the Quasar Accretion Disk Size – Black Hole Mass Relation, in which we confirm that the accretion disk size predictions from simple thin disk theory are too small.

## 1 Introduction

Gravitationally lensed quasars provide a wealth of resources to observers. Their utility in cosmology was realized quite early on (e.g. Refsdal, 1964), and
a number of collaborations (e.g. COSMOGRAIL, Courbin et al. (2005) & H0LiCOW^{1}^{1}1http://www.h0licow.org) continue to pursue measurements of lensed quasar time delays
to make independent measurements of the Hubble Constant, (e.g. Kochanek, 2002; Vuissoz et al., 2008; Fohlmeister et al., 2013; Suyu et al., 2017; Bonvin et al., 2017)
and a range of other useful cosmographic measurements. Quasar microlensing,
also predicted quite some time ago (e.g. Chang & Refsdal, 1979), provides additional motivation for monitoring lensed quasars since the analysis of microlensing variability (e.g. Kochanek, 2004; Morgan et al., 2006; Poindexter et al., 2007; Morgan et al., 2010; Hainline et al., 2013; MacLeod et al., 2015) and chromatic flux ratio anomalies (e.g. Pooley et al., 2007; Bate et al., 2008, 2011; Mediavilla et al., 2011; Pooley et al., 2012; Schechter et al., 2014)
can be analyzed to probe the central engines of the quasars and the properties of the lens galaxy.

The quadruply lensed quasar WFI J2033-4723 (hereafter WFI2033; 203342.08, -4723’43.0” [J2000.0]) was discovered during a wide-field imaging survey for lensed quasars in the southern hemisphere using the MPG/ESO 2.2 m telescope (Morgan et al., 2004). It has a source redshift of = 1.66, a lens redshift of (Eigenbrod et al., 2006) and a maximum image separation of . Vuissoz et al. (2008, hereafter V08) used three seasons (2004-2007) of monitoring data from the Small and Moderate Aperture Research Telescope System (SMARTS) 1.3m telescope at the Cerro Tololo Inter-American Observatory (CTIO) and the 1.2m Leonhard EULER Swiss telescope at La Silla, Chile to measure time delays of days and days between the merged A1+A2 = A, B and C images. V08 found no evidence of variability due to extrinsic factors such as microlensing. Recently, however, Giannini et al. (2017) made a robust detection of microlensing in their 4-season monitoring campaign using the 1.54m Danish telescope at La Silla, a result which we independently corroborated in this investigation. Most recently, Motta et al. (2017) used the single-epoch chromatic microlensing technique to make estimates of the size of the central engine and broad line region in WFI2033.

In this paper, we combine 9 new seasons of WFI2033 monitoring data with the 4 seasons of data from V08 to create a 13-season set of light curves. We present our observational data and reduced light curves in §2, and we analyze the full combined light curves in §3 to confirm the V08 time delays and measure the A1–A2 delay for the first time. In §4 we describe our microlensing analysis technique to include the properties of our strong lens models for WFI2033. In §5, we present the results of our analysis, and we discuss their implications for accretion disk theory. Throughout our discussion we assume a flat cosmology with , , and (Hinshaw et al., 2009).

## 2 Observational Data

### 2.1 Imagery

We observed WFI2033 in the - (F555W), - (F814W) and - (F160W) bands using the Hubble Space Telescope (HST) as an element of the the CfA-Arizona
Space Telescope Survey (CASTLES^{2}^{2}2http://cfa.harvard.edu/castles/, Lehár et al., 2000). The - and -band images were taken using the
Wide-Field Planetary Camera 2 (WFPC2).
The band images, originally presented in V08, were taken using the Near-Infrared Camera and Multi-Object
Spectrograph (NICMOS). We fit the astrometry and photometry of the lens in the HST imagery
with the imfitfits (Lehár et al., 2000) routine, using a de Vaucouleurs model for the lens galaxy G1,
an exponential disk model for the quasar host galaxy and point sources for the quasar images.
Our astrometric and photometric fits, consistent with those made independently by V08, are presented in Table 1.

Component | Astrometry | Photometry | |||
---|---|---|---|---|---|

H=F160W | I=F814W | V=F555W | |||

A1 | |||||

A2 | |||||

B | |||||

C | |||||

G1 |

### 2.2 Monitoring Observations

On the 1.3 m SMARTS telescope we used the optical channel of the dual-beam ANDICAM instrument (DePoy et al., 2003),which has a plate scale of 0369 pixel and a field of view. The mean sampling of the SMARTS data is one epoch every eight days, with three 300 s exposures at each epoch using the R-band filter. The R-band filter has an effective wavelength of 658 nm, translating to a rest-frame wavelength of 2473 Å in the UV. The SMARTS dataset consists of 117 epochs between March 2004 and August 2017.

On the 1.2 m EULER telescope we used the EulerCAM camera which has a plate scale of 02149 pixel and a field of view. The mean sampling of the EULER data is one epoch every five days, with five 360s exposures at each epoch using the or ‘Rouge Genève’ filter. The filter is a modified broad Gunn R filter, with an effective wavelength of 660 nm translating to a rest-frame wavelength of 2481 Å. The new EULER dataset consists of 178 epochs between October 2010 and December 2016.

The details of our photometric measurement technique are
discussed in Kochanek et al. (2006), but we provide a brief summary of that process here. We use five reference stars, located at ,
, , and , with respect to image A1. These reference stars are
used as the basis for a
three-component elliptical Gaussian point-spread function (PSF) model, which we apply to the blended quasar images to obtain the relative brightness of each component
at each epoch.
When applying the PSF model, we hold the relative positions of the lens components fixed to the astrometry from our *HST* -band images.
We model the lens galaxy using a nested Gaussian with fixed effective radius and flux to approximate a de Vaucouleurs profile. For the effective radius, we used
our measurement from the fits , and for the flux we use the
value which minimizes the total in the residuals following galaxy model subtraction when summed over all epochs.
We measured a very small color offset of 0.002 magnitudes between the EULER and SMARTS photometry which we applied to the
SMARTS data when creating our combined light curves and the data provided in Tables 3 and 4.
Both the EULER and SMARTS images
are characterized by a median stellar FWHM (seeing) of . Since the merging A1/A2 pair are separated by only , deconvolving the
flux from these two images was challenging. For seeing conditions worse than and for SMARTS and EULER, respectively, we were unable to reliably resolve any of the
quasar’s images, so we were forced to discard images taken under these conditions. We also discarded 31 of the 326 total observing
epochs from SMARTS and EULER due to bright sky or
cloudy observing conditions.
In Figure 2 we display our new light curves alongside the published light curves from V08. Since V08 were unable to reliably separate the
flux from images A1 and A2, they summed the flux from this closely spaced merging pair to create a single image A light curve in which .

## 3 Time Delays

Using the polynomial light curve fitting technique of Kochanek et al. (2006), we measured the time delays between the combined image A = A1 + A2, image B and image C. In the Kochanek et al. (2006) method, the intrinsic and extrinsic variability in the light curves are fitted by Legendre polynomials, and the polynomial order is chosen using the F-test. In the case of the delay between images A and B and images C and B, we found that a order polynomial provided a sufficient fit for the source variability and that a order polynomial was appropriate for approximating and removing the microlensing variability. In Figure 3, we show the statistic for the time delay fits. The delay measurements and their uncertainties are days and days (in the sense that image B leads both images A and C). In Figure 4, we show the light curves shifted by these delay values. These new measurements are fully consistent with those of V08.

Using our newly reduced light curves, we also obtain a tentative measurement of the delay between images A1 and A2. With a and order polynomials for the intrinsic and microlensing variability, respectively, we find that image A1 leads image A2 by days. We display these results in Figure 5. While V08 were not able to measure the A1-A2 delay, they did constrain the expected range of that delay to days using a series of lens galaxy mass models. While significantly coarser than our measurement of the B-A and B-C delays, the A1-A2 measurement is consistent with the V08 lens models, although this pair will have the largest fractional uncertainties from microlensing-induced variability (Tie & Kochanek, 2018). In the present paper, we generate a series of lens galaxy models in which the expected A1-A2 delay is days, also consistent with our A1–A2 measurement. With these updated time delay measurements, we proceed to the primary goal of this investigation, the analysis of extrinsic variability from microlensing in the reduced light curves. A full analysis of the updated delays will be published in Bonvin et al. (2018).

## 4 Microlensing Analysis

### 4.1 Lens Galaxy Models and Magnification Patterns

In essence, our Bayesian Monte Carlo technique for microlensing analysis is an attempt to reproduce the observed microlensing variability using a large set of models for the physical conditions that might have led to this variability (Kochanek, 2004). All of this hinges on our ability to accurately model the conditions in the lens galaxy through which the quasar’s light must pass. We started by applying the LENSMODEL code of Keeton (2001) to the astrometry from our observations to yield a range of models for the stellar and dark matter content in the lens galaxy at the positions of the lensed images. Following V08, we adopted a 2-component model for the lens galaxy (G1 in Fig. 1). Since this system is now known to exhibit microlensing of both the continuum and the Broad Line Region (BLR) (Sluse et al., 2012; Motta et al., 2017), we did not use the flux ratios or those from our ground-based observations as a constraint on the lens galaxy mass models. We required our models to reproduce the astrometry of the lensed images, and we allowed the position, effective radius, ellipticity and position angle of the lens galaxy to vary within the uncertainties of the photometric model presented in §2.1. Consistent with V08, we were unable to model the astrometry of the lensed images unless we included the influence of the neighboring galaxy G2, the east-west shear from which cannot be created by G1 since it has an ellipticity position angle of only east of north. We modeled G2 as a singular isothermal sphere whose properties were also allowed to vary within the uncertainties of our photometric and astrometric fits.

Since the dark matter content in the lens galaxy is unknown, we created a series of 10 models for the lens galaxy in which the dark matter fraction varies across an order of magnitude. We began by modeling the lens galaxy using only a de Vaucouleurs profile. In each subsequent model, we decreased the monopole moment of the stellar de Vaucouleurs component by 10% of the constant model mass, and we added an extended, concentric Navarro, Frenk, & White (1996, NFW) component to model the dark matter. We parameterized this series using the quantity , representing the fraction of the lens galaxy mass relative to the constant mass-to-light (M/L) ratio model. From this model sequence, we extract the total convergence , the convergence from the stars , the shear and the shear position angle at the location of each lensed image. While models in the range are more consistent with our measured time delays, for completeness we use the entire model sequence in our Monte Carlo microlensing simulations because Schechter & Wambsganss (2002) demonstrated that local microlensing statistics are very sensitive to . The parameters of all 10 models are presented in Table 2.

Convergence | Shear | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

A1 | A2 | B | C | A1 | A2 | B | C | A1 | A2 | B | C | ||

Note. – Convergence , shear and the fraction of the total surface density composed of stars at each image location for the series of mass models. The parameter is the mass of the de Vaucouleurs model for the visible lens galaxy relative to its mass in the absence of dark matter. The per degree of freedom for each model is provided in the column.

Using the inverse ray-shooting technique (as described in Kochanek et al., 2006), we generated forty random realizations of the expected microlensing magnification conditions in the vicinity of each image for each of our 10 macro models parameterized by . The magnification patterns are pixels representing a projected source plane scale of twenty Einstein radii or cm. This implies a pixel scale of cm in the source plane. We assumed an initial stellar mass function (IMF) of with a dynamic range of 50, which approximates the microlensing-based Galactic bulge IMF of Gould (2000), although Wyithe et al. (2000) and Congdon et al. (2007) show that microlensing statistics are not especially sensitive to choice of IMF.

### 4.2 Monte Carlo Method

Armed with 400 sets of magnification patterns for a range of lens galaxy mass models, we used the Monte Carlo light curve fitting technique of Kochanek (2004) as modified by Morgan et al. (2008) to model unresolved image pairs. Since the light curves are a full 13 seasons in length, we binned them using a window of days to reduce computation time. The twenty day binning window was sufficiently short to avoid overly smoothing the microlensing variability while adequately reducing the run time for the Monte Carlo routine. The date of each twenty-day bin was set as the mean Heliocentric Julian Date (HJD) of the measurements included in that bin. Since the light curves from V08 do not provide individual measurements for the A1 and A2 images, we adjusted the statistical weight of the A = A1+A2 data points to appropriately account for the combined fluxes in these cases. Also, for the combined cases we used either the A1 or A2 magnification pattern (with equal frequency) when creating a simulated light curve for the combined A1 + A2 image.

Prior to each Monte Carlo trial, we convolved each set of magnification patterns with a Gaussian surface brightness profile at a range of trial source sizes

(1) |

where is the radius of the accretion disk scaled by , the mean mass of a lens galaxy star. Although we used a Gaussian profile, the exact choice of photometric emission model is unimportant, since Mortonson et al. (2005) showed that microlensing statistics are largely a function of the half-light radius of the emitting region, not the exact properties of the emission profile. In a given trial, a convolved pattern is run past a model point source on a random trajectory and at a random transverse speed 10 km/s km/s. Changes in magnification with time are logged at the epochs of the observed data and a running comparison with the light curves is made. The quality of the fit is tallied in real time using a statistic, and, to save computational time, fits with are discarded since they will not contribute significantly to the Bayesian integrals in our post-run analysis. During the curve fitting process, we allowed for 0.07 and 0.03 magnitudes of systematic error in the photometry of images A1 & A2 and B & C, respectively. We also allowed for 0.5 magnitudes of uncertainty in the intrinsic flux ratios between the lensed images, since both the continuum and the broad line region are affected by microlensing in this lens (Motta et al., 2017) and to allow for the influence of substructure in the lens. We attempted fits per set of magnification patterns for a grand total of trials, requiring approximately two weeks of run time on the US Naval Academy High Performance Computing (HPC) cluster. In Figure 6, we display two of the best fits from our Monte Carlo analysis to the time-delay corrected difference light curves of WFI2033. Consistent with the findings of Giannini et al. (2017), we easily see mags of microlensing variability in the difference light curves from image C. The microlensing in images A1, A2 and B is less pronounced, with mags of extrinsic variability over the 13 seasons of monitoring.

### 4.3 Bayesian Analysis of Monte Carlo Results

Using Bayes’ theorem, the probability of the parameters given the data is

(2) |

where is the collection of physical variables, is the collection of trajectory variables and and are the prior probabilities for the trajectory and physical variables, respectively. The physical variables are parameters of the local magnification tensors (mean convergence and mean shear ), the local properties of the stars (surface density of stars , mass of the average microlens ), the scale radius of the source , and the effective velocity of the source . The probability distribution for any variable of interest can be obtained by marginalizing over the other variables of the simulation. For example, to find the probability density for the scale radius ,

(3) |

where is the probability of fitting the data in a particular trial, sets the priors on the microlensing variables & , and is the (uniform) prior on the scale radius. The total probability is then normalized so that .

We initially do the calculation in “Einstein units”, where all lengths depend upon a factor of the unknown mean stellar mass, . For example, in Figure 8, we display the probability density for the accretion disk scale radius in Einstein Units in which the plotted values assume . This degeneracy can be broken, however, by examining Figure 7, where we display the probability density for the scaled effective velocity (Einstein Units) from the Monte Carlo simulation. We also display a model in physical units for the expected transverse velocity which serves as the statistical prior on . Since , we convolve the prior on with the probability density for to produce a probability density for . The probability density for is then used to convert all scaled lengths (e.g. ) into true, physical units (e.g. ) by convolving with the quantity of interest.

We construct the prior on the transverse velocity, , using the method described in Kochanek (2004). For the peculiar velocity components of both lens and source, we make redshift-based estimates from the models of Mosquera & Kochanek (2011). We estimate the velocity dispersion of the lens galaxy from its Einstein radius, assuming the galaxy is a singular isothermal sphere with relaxed dynamics, which Treu & Koopmans (2004) and Bolton et al. (2008) show is a good approximation.

We display the probability density for the scale radius in physical units in Figure 8. In this plot, we also show a probability density for obtained by assuming a uniform prior on the median lens galaxy stellar mass of . A brief inspection of the plot reveals that the results without a prior on the microlens mass are robustly consistent with the results using the uniform mass prior. As a final step we must correct the scale radius for the disk’s inclination by multiplying by , which is necessary because we have assumed a face-on disk in our simulations and microlensing amplitudes depend on the projected area of a source rather than the shape. We adopt the measurement made without the mass prior at , where is the inclination angle.

## 5 Results & Discussion

In Figure 9, we plot the size of the accretion disk in WFI2033 on the Accretion Disk Size - Black Hole Mass Relation (Morgan et al., 2010) assuming an inclination angle , where we have corrected the scale radius at the wavelength corresponding to the center of the rest-frame -band, Å, to , the scale radius at Å, assuming the scaling of Shakura & Sunyaev (1973) thin disk theory. We assume or , the expectation value for the inclination of a randomly oriented disk. For the black hole mass (), we use the result from Sluse et al. (2012), who used the Mg ii emission line to find . Motta et al. (2017) also estimated in this system using Keplerian dynamics, but their method yielded very large uncertainties with . The Morgan et al. (2010) relation was derived using estimates based on emission line widths so the emission line width-based measurement from Sluse et al. (2012) is a better choice.

The expectation value for the scale radius at without the prior on the mass of microlenses is . This is fully consistent with the results of Motta et al. (2017), who estimate a scale radius cm at Å using single-epoch spectroscopy, which, when scaled to 2481 Å assuming is . The Motta et al. (2017) result is strongly dependent upon priors, especially the assumption of a median microlens mass , nevertheless, the independence of the techniques provides robust support for our result. Blackburne et al. (2011) also estimated the scale radius of the accretion disk in WFI2033 using single epoch, multi-wavelength photometry, but their results are, by their own admission, anomalous, as they also predict with highest likelihood an accretion disk with an inverted (increasing toward the outer edge) temperature profile. Like Motta et al. (2017), our size measurement is significantly smaller than the Blackburne et al. (2011) estimate.

In Figure 9, we made several updates in addition to the new WFI2033 measurement. Measurements of the accretion disk scale radius using the microlensing variability technique of Q 0957+561 (Hainline et al., 2012), SBS 0909+532 (Hainline et al., 2013) were added, and updates to the QJ 0158–4723 (Morgan et al., 2012), and SDSS 0924+0219 (MacLeod et al., 2015) measurements were also included. With changes to 2 out of the 11 existing points and the addition of 3 new measurements, the updated Accretion Disk Size - Black Hole Mass Relation (Morgan et al., 2010) is

(4) |

This is consistent with the original fit from Morgan et al. (2010), cmM, and the shallower slope brings the relation into excellent agreement with the expectation from thin disk theory ().

There are now 14 systems in which the accretion disk size has been measured using the microlensing variability (e.g. Kochanek, 2004) technique. With the exception of SBS 0909+532, in which the luminosity-based size estimate is marginally larger than the microlensing-based size measurement, microlensing-based size measurements are consistently larger than the luminosity-based thin disk size estimates by an average of dex. The SBS 0909+532 luminosity-based thin disk size estimate is somewhat suspect, however, since Sluse et al. (2012) and Lehár et al. (2000) found very different photometric fits for the lens galaxy in this system, leading to significant uncertainty in the magnification and, consequently, the intrinsic luminosity. Very recent continuum emission region reverberation mapping studies in local, lower luminosity AGN (e.g. Fasnaugh et al., 2018; Edelson et al., 2017; McHardy et al., 2016) have revealed similar discrepancies between observed accretion disk size measurements and the predictions of thin disk theory. In Morgan et al. (2010), we proposed that real accretion disks lack the necessary surface brightness to produce their observed luminosity from the smaller area of a simple thin disk model, and we remain confident in that conclusion. We were nevertheless intrigued to find that the slope of vs is remarkably consistent with the predictions of thin disk theory (), so it is the intercept in the accretion disk-size black hole mass relation that is inconsistent with thin disk theory, rather than the slope. In Morgan et al. (2010), we suggested that the most promising explanation for the discrepancy is that accretion disks may have shallower temperature slopes than that predicted by thin disk theory , and this hypothesis remains fully viable. We are hopeful that our ongoing lensed quasar monitoring campaign in the infrared (, and band), corresponding to optical emission in the rest frame, will allow for measurements of accretion disk temperature profiles in the near future.

HJD | A1 | A2 | B | C | Source | |
---|---|---|---|---|---|---|

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS | ||||||

SMARTS |

Note. – HJD is the Heliocentric Julian Day –2450000 days. The goodness of fit of the image, , is used to rescale the formal uncertainties by a factor of . The Image A1-C columns give the magnitudes of the quasar images relative to the comparison stars.

HJD | A1 | A2 | B | C | Source | |
---|---|---|---|---|---|---|

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||

EULER | ||||||