text:
# CHANDRA OBSERVATIONS OF THE CLOVERLEAF QUASAR H 1413+117: A UNIQUE LABORATORY FOR MICROLENSING STUDIES OF A LoBAL QUASAR

## 1 Introduction

Recently, the gravitationally lensed (GL) quasar H 1413+117 (the Cloverleaf quasar) was observed for $``$ 40 ks with the Chandra X-ray Observatory (Oshima et al. 2001, Gallagher et al. 2002). The original goal of this observation (PI: M. Hattori) was to detect a cluster of galaxies that is thought to contribute to the lensing of this system (Kneib et al. 1998; Chae & Turnshek 1999). X-ray emission from the cluster lens was not detected. However, a remarkable iron emission feature at $`E`$=6.21$`\pm `$0.16 keV (90% confidence) was discovered in the combined spectrum of all lensed images. This redshifted emission feature was interpreted as Fe K$`\alpha `$ fluorescence from the far side of the quasar outflow.

H 1413+117 is a gravitationally lensed system consisting of four images of a distant ($`z`$=2.55) broad absorption line quasar (Hazard et al. 1984; Magain et al. 1988). $`\mathrm{𝐻𝑆𝑇}`$ NICMOS-2 observations indicate that the lens is located approximately at the geometrical center of the four images. A firm spectroscopic redshift of the lens has yet to be obtained, however, a cluster of galaxies at a redshift of $`z`$ = 1.7 has been reported to contribute to the lensing of this system (Kneib et al. 1998). H 1413+117 was the first quasar to be detected in the submillimeter continuum and in carbon monoxide emission (Barvainis et al. 1995). H 1413+117 is also a low-ionization broad absorption line quasar (LoBAL) showing low-ionization lines of Al III and Fe III (Monier et al. 1998). LoBALs are BAL quasars that show one or more of the low-ionization lines of Mg II, Al III, Fe II and Fe III (e.g., Boroson & Meyers 1992). Its first X-ray detection was made with the ROSAT PSPC (Chartas et al. 2000). X-ray observations of BAL quasars indicate that they are heavily absorbed with measured column densities ranging between 1021 and 1024 cm-2. The extreme X-ray weakness of LoBALs, their low values of $`\alpha _{\mathrm{OX}}`$, and their high levels of UV continuum polarization are taken to imply that the absorbing gas in the direct line of sight to the quasar central engine is Compton-thick to X-ray radiation from the quasar nucleus and only scattered or reflected radiation can reach the observer along an indirect line of sight (Green et al. 2001; Gallagher et al. 2002).

Here, we present a re-analysis of the *Chandra* observation of H 1413+117 that places constraints on the spatial structure of the X-ray iron line and continuum scattering regions of this LoBAL quasar. The main difference between our and previous analyses is that we performed a spatial and spectral analysis of the individual X-ray images whereas previous X-ray studies were based on the combined spectrum of all images. Throughout this paper we adopt a cosmology in which $`H_0`$ = 75 $`h_{75}`$ km s-1 Mpc-1 and $`q_0`$ = 0.5.

## 2 Observations and Data Analysis

H 1413+117 was observed with *Chandra* for about 38 ks on 2000 April 19. Details of this observation were described in Oshima et al. (2000) and Gallagher et al. (2002). For preparing the *Chandra* data for analysis we used the CIAO 2.2 and CALDB 2.12 products provided by the *Chandra* X-ray Center (CXC). To improve the spatial resolution we removed a $`\pm `$ 0$`^{\prime \prime }.`$25 randomization applied to the event positions in the CXC processing and employed a sub-pixel resolution technique developed by Tsunemi et al. (2001) and later improved by Mori et al. (2001). The deconvolved Chandra ACIS image of H 1413+117 is shown in Figure 1. A comparison between the deconvolved X-ray image and the optical HST images of H 1413+117 taken from the CfA-Arizona Space Telescope Lens Survey (CASTLES)111HST images of H 1413+117 in the F160W, F702W and F814W HST bands are available from the Castles web site: http://cfa-www.harvard.edu/glensdata/Individual/H1413.html. clearly shows the enhancement of X-ray image A compared to that of the optical image A. The X-ray image positions are in agreement within errors with those of the HST images indicating that the X-ray emission used in our spatial-spectral analysis originates from the lensed quasar and not from extended emission of the lensing galaxy or host galaxy. To account for the recently discovered decay in the quantum efficiency of ACIS (possibly caused by molecular contamination of the ACIS filters), we have applied a time-dependent correction implemented in the `XSPEC` model `ACISABS1.1`.222ACISABS is an XSPEC model contributed to the *Chandra* users software exchange world wide web-site http://asc.harvard.edu/cgi-gen/cont-soft/soft-list.cgi. The ACIS quantum efficiency decay is insignificant at energies above 1 keV and does not affect the main results of our analysis. To estimate the X-ray flux ratios we modeled the *Chandra* images of A, B, C and D with point spread functions (PSF’s) generated by the simulation tool `MARX` (Wise et al 1997). The X-ray event locations were binned with a bin-size of 0$`^{\prime \prime }.`$0246. The simulated PSFs were fitted to the *Chandra* data by minimizing the Cash $`C`$-statistic formed between the observed and simulated images of H 1413+117. The relative positions of the images were fixed to the observed NICMOS values taken from CASTLES. We find that the X-ray flux ratios in the full (0.2-8keV) band are \[B/A\]full = 0.49 $`\pm `$ 0.08, \[C/A\]full = 0.35$`\pm `$ 0.07, \[D/A\]full = 0.37 $`\pm `$ 0.07 and the number of detected X-ray events in images A, B, C, and D were 147 $`\pm `$ 13, 72 $`\pm `$ 10, 52 $`\pm `$ 8, and 54 $`\pm `$ 9, respectively. These flux ratios are consistent with the values obtained from the deconvolved X-ray image shown in Figure 1. For comparison the HST R-band flux ratios are \[B/A\]R = 0.94 $`\pm `$ 0.01, \[C/A\]R = 0.78 $`\pm `$ 0.01, \[D/A\]R = 0.74 $`\pm `$ 0.01. The R-band magnitudes were taken from the CASTLES Web site. The flux fraction of image A, \[A/(A+B+C+D)\], is a factor of 1.55 $`\pm `$ 0.17 larger in the full X-ray band than in the HST R-band.

The earlier spectral analyses of H 1413+117 used the combined spectrum of all images. In our re-analysis of the data we search for possible spectral variations between images by extracting a spectrum for image A and a combined spectrum of images B, C, and D. The low S/N of these data did not allow us to extract individual spectra for images B, C, and D. The spectrum of image A was extracted from a circular region of radius 1′′centered 0$`^{\prime \prime }.`$5 SW from the center of image A to minimize contamination from the other images. The combined spectrum of images B, C, and D was extracted from a circular region centered on H 1413+117 with a radius of 2$`^{\prime \prime }.`$5, excluding the source region of image A. We determined corrections to the ancillary response files by simulating the spectra of point sources at the locations of the images within the apertures used in our analysis. For our simulations we used `XSPEC` to generate the source spectra and the ray-tracing tool `MARX` to model the dependence of photon scattering with energy. The background was extracted from a source-free annulus centered on H 1413+117 with inner and outer radii of 10′′ and 30′′, respectively. All errors reported below are at the 90% confidence level unless noted otherwise. A variety of models were fitted to the spectrum of image A and the combined spectrum of images B, C, and D, employing the software tool `XSPEC v11.2` (Arnaud 1996). Spectral fits were performed in the 1–8 keV bandpass. Given the relatively low-S/N of the spectra we performed fits using the $`C`$-statistic (Cash 1979) and compared the results with those obtained using the $`\chi ^2`$ statistic. Spectral fits performed with the $`C`$-statistic do not require binning of the data and produce reliable estimates of the best-fit parameters and their errors for low-S/N spectra. One of the limitations of the $`C`$-statistic is that it cannot provide a goodness-of-fit measure. However, simple analytic expressions do exist for estimating the goodness-of-fit when using the $`\chi ^2`$ statistic. We therefore used fits performed with the $`\chi ^2`$ statistic to evaluate the probabilities of exceeding the derived $`\chi ^2`$ values. For the spectral fits performed with the $`\chi ^2`$ statistic we binned the data to have a minimum of twelve counts per bin. We did not find any significant differences between the estimated best fit parameters and their associated errors provided by the $`C`$-statistic and the $`\chi ^2`$ statistic. In Table 1 we list the best fit parameters and their errors based on fits using the $`C`$-statistic.

The spectra were initially fitted individually with a model consisting of a simple power law with Galactic absorption due to neutral gas with a column density of $`N_\mathrm{H}`$ = 1.8 $`\times `$ 1020 cm-2 (Stark et al. 1992). The model also included neutral intrinsic absorption at $`z=2.55`$. The fit residuals in image A showed an emission feature between 1.65–1.85 keV that led to an unacceptable fit, with a $`\chi ^2`$ = 13.56 for 10 degrees of freedom (dof). Only the spectrum of image A shows clear evidence of this strong emission line. We modeled the residual feature in image A by considering a Gaussian component near 1.7 keV in the observed frame. This led to a significant improvement in fit quality at the $``$ 98% confidence level (according to the $`F`$-test, $`F=5.6`$) with a $`\chi ^2`$ = 4.0 for 7 dof. In Figure 2a we show the spectrum of image A. Protassov et al. (2002), have recently argued that the $`F`$-test cannot be applied when the null values of the additional parameters fall on the boundary of the allowable parameter space. They proposed a Monte Carlo approach to determine the distribution of the $`F`$ statistic. We followed this approach and constructed the probability density distribution of the $`F`$-statistic between spectral fits 1 and 2 of Table 1. Our Monte Carlo simulations indicate that the probability of obtaining an $`F`$ value of 5.6 or greater is $``$1.5 $`\times `$ 10-2 close to the result found by using the analytic expression for the distribution of the $`F`$-statistic. The best-fit rest-frame energy of the emission feature is $`E_{line}`$ = 6.3$`\pm 0.1`$ keV. Adding an unresolved Gaussian line to the combined spectrum of images B, C, and D leads to a slight improvement in fit quality at the 86% confidence level (according to the $`F`$-test). The combined observed-frame spectrum of images B, C, and D fit with Galactic absorption, neutral absorption at the source, and a power-law model is shown in Figure 2b. The best-fit values for the photon indices in the spectrum of image A and the combined spectrum of images B, C, and D, are $`\mathrm{\Gamma }`$ = 1.8$`\pm `$0.8 and $`\mathrm{\Gamma }=1.7\pm 0.7`$, respectively (see fits 1 and 4 of Table 1). Our fits support the presence of an intrinsic absorber towards image A with a column density of $`N_\mathrm{H}`$ = (2.3$`\pm `$2.2) $`\times `$ 1023 cm-2 (see fit 1 of Table 1). We also find a significant intrinsic absorber from the fit to the combined spectrum of images B, C, and D of $`N_\mathrm{H}`$ = (3$`{}_{}{}^{+3}{}_{2}{}^{}`$) $`\times `$ 1023 cm-2. The observed column in images B, C, and D (see fit 4 of Table 1) may represent the column density perpendicular to the outflow rather that along the outflow as shown in Figure 3 of Oshima et al. (2001). They estimated the column along the flow to be $`N_H10^{24}`$ cm-2. More complex models including ionized absorption or partial covering were not considered in our analysis due to the low S/N of the available spectrum and the high $`z`$ of H 1413+117 which results in several absorption edges due to Ne, O, N, and C being redshifted outside the *Chandra* bandpass.

We performed simultaneous fits to the spectrum of image A and the combined spectrum of images B, C, and D to estimate the significance of the difference in the Fe K$`\alpha `$ lines. The common spectral model for the fit consisted of (a) Galactic absorption due to neutral gas with a column density of $`N_\mathrm{H}`$ = 1.8 $`\times `$ 1020 cm-2, (b) neutral absorption at the redshift of the quasar, (c) a simple power law, and (d) a Gaussian component near 1.7 keV in the observed frame. The photon index $`\mathrm{\Gamma }`$ and the energies and widths of the Fe K$`\alpha `$ lines were linked as common parameters in this spectral fit. The best-fit parameters and 90% confidence errors are presented in Table 1. The rest-frame equivalent widths and 90% confidence levels of the lines in the spectrum of image A and the combined spectrum of images B, C, and D are $`EW_A`$ = 1.27$`{}_{}{}^{+0.59}{}_{0.59}{}^{}`$ keV and $`EW_{BCD}`$ = 0.20$`{}_{}{}^{+0.36}{}_{0.20}{}^{}`$ keV, respectively. In Figure 3 we show the 68.3% and 85% confidence contours of the normalizations of the Fe K$`\alpha `$ lines versus the rest-frame Fe K$`\alpha `$ line energy. We chose to show the 85% contours instead of the commonly used 90% contours because the 85% contours touch exactly in this case. Based on these contours, the probability that the iron lines in the spectrum of image A and the combined spectrum of images B, C, and D have strengths that are consistent with each other is less than 0.02. For this confidence contour plot we fixed the parameters of the power-law component of the model to their best-fit values listed in fit 5 of Table 1. We conclude that the Fe K$`\alpha `$ line in the spectrum of image A is significantly enhanced by a factor of about six from that of the combined spectrum of images B, C, and D.

## 3 Discussion

Our re-analysis of the X-ray spectra of the individual lensed images of H 1413+117 indicates the presence of a microlensing event in image A because (a) the flux fraction of image A, \[A/(A+B+C+D)\], is a factor of 1.55 $`\pm `$ 0.17 larger in the X-ray band than what was detected in the HST R-band, (b) a redshifted narrow iron line is present with a relatively high significance only in the spectrum of image A. To obtain additional insight on the possible energy dependence of the magnification of image A, we estimated relative enhancements of image fluxes in the continuum and iron-line emission components based on fit 5 of Table 1. Specifically, we find that the ratio of the X-ray fluxes of image A and the combined images B, C and D in the iron line component is \[A/(B+C+D)\]$`_{\mathrm{Fe}\mathrm{K}\alpha }`$ $`>`$ 1.7 at the 85% confidence level (see Figure 2). The ratio of the X-ray flux densities at 1 keV (observed frame) in the continuum component is \[A/(B+C+D)\]Cont = 0.8 $`\pm `$ 0.2 (90% confidence intervals). These values differ significantly from the HST R-band flux ratio of \[A/(B+C+D)\]R=0.41$`\pm `$0.01 which is expected to be less affected by a microlensing event produced by a star or group of stars in the foreground lensing galaxy since the characteristic sizes of the X-ray continuum emission regions of quasars are expected to be considerably less than the sizes of the optical R-band emission regions. Specifically, variability studies indicate that the sizes of the X-ray and optical continuum emission regions are $``$ 4$`\times 10^{1516}`$ cm and $``$ 4$`\times 10^{1617}`$ cm, respectively (e.g., Chartas et al. 2001; Dai et al 2002; Wyithe et al. 2000). We note that since the flux fraction of image A in the X-ray continuum is greater than the flux fraction in the R band, the large equivalent width of the Fe line in image A is not caused by a decrease in the X-ray continuum but is more likely produced by an amplification of the flux of the line.

A scenario in which intrinsic variability is the cause of the detection of a strong iron line in only image A and a significantly weaker line in the combined spectrum of images B, C and D requires fine tuning such that the flare of the iron line needs to occur within a certain window in time for it to be detected with a large $`EW`$ in only one image and the flare must last for only a few days. Our simple two singular isothermal sphere lens model for H 1413+117 predicts that image B leads images C, A, and D by 5 days, 9 days and 16 days respectively. A similar lens model by Chae & Turnshek (1999) was found to reproduce the observed relative amplifications fairly accurately. Another strong argument against the intrinsic variability interpretation is that changes of factors of $``$ 6 in the equivalent width of a fluorescent Fe K$`\alpha `$ line in a Seyfert 1 galaxy over periods of a few days have not been detected to date (e.g., Markowitz, Edelson & Vaughan 2003; Gliozzi, Sambruna & Eracleous 2003; Lee et al. 2000; Leighly 1999). We anticipate that the amplitude of Fe K$`\alpha `$ line variability to be less in quasars than in Seyfert 1 galaxies due to the larger sizes and luminosities of quasars. X-ray variability, at least of the continuum emission in Seyfert 1 galaxies and quasars, is found to be anticorrelated with X-ray luminosity (e.g., Nandra et al. 1997; Turner et al. 1999; George et al. 2000; Dai et al. 2003). A more likely explanation of the observed energy dependent magnification in the X-ray band and the enhancement of image A in the X-ray band compared to that in the optical is a microlensing event in image A. Such an event will not lead to time-delayed magnifications in the remaining images and therefore can explain the non-detection of the Fe K$`\alpha `$ line in the remaining images. A characteristic size scale commonly used in microlensing analysis is the projected Einstein ring radius of the star, $`\zeta _E`$. Emission regions with a size significantly larger than $`\zeta _E`$ will not be affected by microlensing, whereas, emission regions less than $`\zeta _E`$ will be significantly magnified. For H 1413+117 with lens and source redshifts of 1.7 and 2.55 respectively, and assuming an isolated star, the Einstein-ring radius on the source plane is $`\zeta _E`$ $``$ 2 $`\times `$ 1016$`(M_{\mathrm{star}}/M_{})^{1/2}`$$`h_{75}^{0.5}`$ cm.

The observed enhancement of the Fe K$`\alpha `$ line implies that the size of this iron line region is $`R_{line}`$ $`\stackrel{<}{}`$ $`\zeta _E`$. The enhanced Fe K$`\alpha `$ emission observed in image A may be produced by microlensing of Fe K$`\alpha `$ fluorescence from the far side of the outflow. Assuming the Murray et al. (1995) disk-wind geometry for BALQSOs we expect that Fe K$`\alpha `$ fluorescence from the near-neutral accretion disk, thought to be emitted (based on modeling of the Fe K$`\alpha `$ line profiles in several Seyferts) from annuli of inner radii of $`r_{in}`$ $``$$`r_g`$ and outer radii of $`r_{out}`$ $``$ 30 $`r_g`$ (e.g., Nandra et al. 1997), where $`r_g=GM_{bh}/c^2`$ and $`M_{bh}`$ is the mass of the black hole of H 1413+117, will be heavily absorbed by the quasar outflow. The far side of the wind is a more plausible source of the Fe K$`\alpha `$ line because our line of sight to it does not traverse a significant absorbing column. We estimate the size of the microlensed fluorescent region of the quasar outflow to be $`\stackrel{<}{}`$ $`\zeta _E`$. We have assumed a minimum launching radius for the UV BAL absorber of $`r_{\mathrm{min}}1\times 10^{16}(M_{bh}/10^8M_{})^{1/2}`$ cm based on hydrodynamical simulations performed by Proga et al. (2000). To estimate $`M_{bh}`$ we assumed the lensed bolometric luminosity of $`L_{Bol}`$ = 3.9 $`\times `$ 1047 $`\mu ^1`$ erg s-1 (Granato et al. 1996), $`\mu `$ is the magnification factor ranging between 10 and 23 based on lens models of this system (e.g., Venturini & Solomon 2003, Chae et al. 1999), $`L_{Bol}/L_{Edd}`$ = $`\eta `$, $`L_{Edd}`$ is the Eddington luminosity and $`\eta `$ is thought to range between 0.1 and 1 in quasars. We find, $`L_{Edd}`$ = 3.9 $`\times `$ 1047 $`\mu ^1\eta ^1`$ erg s-1, $`r_g`$ = 4.5 $`\times `$ 1014 $`\mu ^1\eta ^1`$ cm, $`M_{bh}`$ = 3 $`\times `$ 109 $`\mu ^1\eta ^1`$ $`M_{}`$ and $`r_{\mathrm{min}}`$= 5.5 $`\times `$ 1016 $`\mu ^{1/2}\eta ^{1/2}`$ cm. We note that with the values we have adopted the product $`\mu \eta `$ is approximately equal to 1. The relatively smaller magnification of the observed X-ray continuum emission compared to the line emission suggests that the size of the region producing the observed X-ray continuum emission is larger than $`\zeta _E`$.

The processes thought to contribute to the continuum X-ray emission in radio-quiet quasars include inverse Compton scattering of photons in disk coronae by UV photons originating in the accretion disk, resulting in an increase in the energy of the upscattered photons from the UV range to the 2-10 keV range. At rest energies above $``$ 10 keV, Compton scattering of photons from the disk corona by the disk becomes significant. Variability studies of radio-quiet quasars with similar black hole masses to H 1413+117 indicate that the size of the emission region responsible for these continuum processes described above is of order 1015-16 cm (e.g., Chartas et al. 2001; Dai et al 2002.) This, combined with our conclusion that the observed X-ray continuum emission must have a size $`>>`$ $`\zeta _E`$ suggests that the X-ray continuum emission in H 1413+117 is not viewed directly but most likely is reflected in a scattering region with a characteristic size of $`r_{scat}`$ $`>>`$ $`\zeta _E`$.

Our constraint on the size of the scattering region is supported by recent HST observations of H 1413+117; using imaging polarimetry with the HST, Chae et al. (2001) detected a 0.07 mag dimming and a change of the linear polarization component of image D relative to the other images during the June 1999 observation which they interpreted as microlensing of a portion of a large-scale, asymmetric continuum scattering region seen in polarized light. They also conclude that the microlensing caustic during their observations was located away from the continuum emitting region. They independently conclude that the scattering region is 10$`{}_{}{}^{18}L_{46}^{0.5}`$ cm $`>`$ $`r_{scat}`$ $`>`$ 2 $`\times `$ 1016$`(M_{\mathrm{star}}/M_{})^{1/2}`$$`h_{75}^{0.5}`$ cm. We also note that based on the simulations of Witt, Mao, and Schechter (1995) images A and D of H 1413+117 have the largest probability of being microlensed due to the relatively large optical depths for microlensing events along these lines of sight. In Figure 4 we present a sketch of a scenario that can explain the observed X-ray spectra of H 1413+117. For simplicity the caustic plane is oriented perpendicular to the observer’s line of sight. In this proposed scenario the microlensing caustic is centered near the black hole and the observed continuum emission in all images originates predominatly from a large scale scattering region. We expect the continuum X-ray emission region and the Fe K$`\alpha `$ fluorescence from the near-neutral accretion disk to be significantly absorbed along our line of sight. This is supported by the fact that H 1413+117 is a LoBAL quasar and recent studies indicate that their X-ray spectra may be reflection/scattered dominated (e.g., Green et al., 2001; Gallagher 2002). Based on the anticorrelation between $`\alpha _{ox}`$ and 2500 $`\mathrm{\AA }`$ luminosity (Vignali et al. 2003) for radio-quiet quasars (RQQs) we estimate an average value of $`\alpha _{ox}`$ $``$ -1.6 $`\pm `$ 0.2 for a RQQ with $`l_{2500\mathrm{\AA }}`$ = 6.1 $`\times `$ 1030 erg s-1 Hz-1 of H 1413+117. In a recent survey of high redshift lensed quasars Dai et al. (2003) calculated $`\alpha _{ox}`$($`\alpha _{ox}`$ corrected) = $`3.3\pm 0.5(1.8_{0.1}^{+0.15})`$ for H 1413+117 based on the X-ray flux of the images B, C, and D (image A was not included since it is likely microlensed) confirming the X-ray weakness of this LoBAL. The estimated values of $`\alpha _{ox}`$ are highly dependent on the assumed model for the intrinsic absorption. The observed X-ray emission in image A also contains significantly magnified Fe K$`\alpha `$ emission originating from the far side of the flow and less magnified emission from a portion of the large scale scattering region.

Another possiblity is that the microlensing caustic during the Chandra observations was not centered on the black hole but located further away than the X-ray continuum emission region of the inner accretion disk. In this scenario the observed continuum emission in all images originates from both the continuum-emitting region near the inner accretion disk and the large-scale scattering region. The X-ray flux observed in image A also contains magnified emission from the far side of the outflow and magnified emission from a portion of the large scale scattering region. Since the microlensing caustic in this case does not overlap with the central region of the black hole complete obscuration of the emission from the inner region of the accretion disk is not required to explain the observed energy-dependent magnification. To summarize, the various emission regions referred to in this analysis in order of smallest to largest regions are: (a) A very small-scale ($`<<\zeta _E`$), possibly obscured Fe K$`\alpha `$ fluorescence region. (b) A small-scale ($`\stackrel{<}{}\zeta _E`$) portion of the reflection region from the far side of the outflow, which is relatively cold and neutral, slightly redshifted, and is highly amplified by microlensing. This region causes the large equivalent width for the iron line in image A. (c) A small-scale ($`\stackrel{<}{}\zeta _E`$) X-ray continuum emission region originating near the inner-accretion disk. (d) A mid-scale ($`\zeta _E`$) optical continuum R band emission region, and, (e) A large-scale ($`>\zeta _E`$) X-ray and optical scattering region which is highly ionized and reflects the X-ray source without creating a fluorescent iron line.

Optical spectra of H 1413+117 taken with the SILFID spectrograph at the Canada-France-Hawaii (CFH) telescope in 1989 (Angonin et al. 1990) and the HST Space Telescope Imaging Spectrograph in 2000 (spectra taken by E. M. Monier cited in Chae et al. 2001) indicate that a microlensing event along the sight-line of image D was in progress for at least a duration of 11 years. This timescale is consistent with the estimate by Witt, Mao, & Schechter (1995) of the time for the source to cross an Einstein radius for H 1413+117. We may, therefore, expect that the possible microlensing event in image A persists for a similar timescale.

The detection of the iron line is possible perhaps because of an ongoing microlensing event in image A. Future observations of H 1413+117 with *Chandra* may provide a confirmation of the microlensing interpretation if a significant variation in the continuum and/or iron line flux ratios in the images is detected. This confirmation would validate our microlensing models and our conclusions regarding the relative size of the continuum and iron line scatterers.

We acknowledge financial support from NASA grants NAG5-9949, NAS 8-38252 and NAS 8-01128. We thank the anonymous referee for helpful comments and suggestions.
## References

> Angonin, M.-C., Vanderriest, C., Remy, M., & Surdej, J. 1990, A&A, 233, L5
>
> Arnaud, K. A. 1996, ASP Conf. Ser. 101: Astronomical Data Analysis Software and Systems V, 5, 17
>
> Barvainis, R., Antonucci, R., Hurt, T., Coleman, P., & Reuter, H.-P. 1995, ApJ, 451, L9
>
> Boroson, T. A. & Meyers, K. A. 1992, ApJ, 397, 442
>
> Cash, W. 1979, ApJ, 228, 939
>
> Chartas G., 2000, ApJ, 531, 81.
>
> Chartas, G., Dai, X., Gallagher, S. C., Garmire, G. P., Bautz, M. W., Schechter, P. L., & Morgan, N. D. 2001, ApJ, 558, 119
>
> Chartas, G., Brandt, W. N., Gallagher, S. C., & Garmire, G. P. 2002, ApJ, 579, 169
>
> Chartas, G., Brandt, W. N., Gallagher, S. C., & Garmire, G. P. 2003, Astron. Nachr., 324, 173
>
> Chae & Turnshek 1999, ApJ, 514 587
>
> Chae, K., Turnshek, D. A., Schulte-Ladbeck, R. E., Rao, S. M., & Lupie, O. L. 2001, ApJ, 561, 653
>
> Dai, X., Chartas, G., Agol, E., Bautz, M. W., & Garmire, G. P. 2003, ApJ, 589, 100
>
> Dai, X., Chartas, G., Eracleous, M., & Garmire, G. P. 2003, submitted to ApJ
>
> Gallagher, S. C., Brandt, W. N., Chartas, G., & Garmire, G. P. 2002, ApJ, 567, 37
>
> George, I. M., Turner, T. J., Yaqoob, T., Netzer, H., Laor, A., Mushotzky, R. F., Nandra, K., & Takahashi, T. 2000, ApJ, 531, 52
>
> Gliozzi, M., Sambruna, R. M., & Eracleous, M. 2003, ApJ, 584, 176
>
> Granato, G. L., Danese, L., & Franceschini, A. 1996, ApJ, 460, L11
>
> Green, P. J., Aldcroft, T. L., Mathur, S., Wilkes, B. J., & Elvis, M. 2001, ApJ, 558, 109
>
> Hasinger, G., Schartel, N., & Komossa, S. 2002, ApJ, 573, L77
>
> Hazard, C., Morton, D. C., Terlevich, R., & McMahon, R. 1984, ApJ, 282, 33
>
> Kneib, J.-P., Alloin, D., & Pello, R. 1998, A&A, 339, L65
>
> Lee, J. C., Fabian, A. C., Reynolds, C. S., Brandt, W. N., & Iwasawa, K. 2000, MNRAS, 318, 857
>
> Leighly, K. M. 1999, ApJ, 125, S317
>
> Magain, P., Surdej, J., Swings, J.-P., Borgeest, U., & Kayser, R. 1988, Nature, 334, 325
>
> Markowitz, A., Edelson, R., & Vaughan, S. 2003, ApJ, in press
>
> Monier, E. M., Turnshek, D. A., & Lupie, O. L. 1998, ApJ, 496, 177
>
> Mori, K., Tsunemi, H., Miyata, E., Baluta, C., Burrows, D. N., Garmire, G. P., & Chartas, G. 2001, in ASP Conf. Ser. 251, New Century of X-Ray Astronomy, ed. H. Inoue & H. Kunieda (San Francisco: ASP), 576
>
> Nandra, K., George, I. M., Mushotzky, R. F., Turner, T. J., & Yaqoob, T. 1997, ApJ, 476, 70
>
> Nandra, K., George, I. M., Mushotzky, R. F., Turner, T. J., & Yaqoob, T. 1997, ApJ, 477, 602
>
> Oshima, T., Mitsuda, K., Fujimoto, R., Iyomoto, N., Futamoto, K., Hattori, M., Ota, N., Mori, K., Ikebe, Y., Miralles, J. M., & Kneib, J.-P. 2001, ApJ, 563, L103
>
> Protassov, R., van Dyk, D. A., Connors, A., Kashyap, V. L., & Siemiginowska, A. 2002, ApJ, 571, 545
>
> Shalyapin, V. N. 2001, Astronomy Letters, 27, 150
>
> Stark, A. A., Gammie, C. F., Wilson, R. W., Bally, J., Linke, R. A., Heiles, C., & Hurwitz, M. 1992, ApJS, 79, 77
>
> Tsunemi, H., Mori, K., Miyata, E., Baluta, C., Burrows, D. N., Garmire, G. P., & Chartas, G. 2001, ApJ, 554, 496
>
> Turner, T. J., George, I. M., Nandra, K., & Turcan, D. 1999, ApJ, 524, 667
>
> Vignali, C., Brandt, W. N., Schneider, D. P., Anderson, S. F., Fan, X., Gunn, J. E., Kaspi, S., Richards, G. T., & Strauss, M. A., 2003, AJ, 125, 2876
>
> Wise, M. W., Davis, J. E., Huenemoerder, Houck, J. C., Dewey, D. Flanagan, K. A., and Baluta, C. 1997, The MARX 3.0 User Guide, CXC Internal Document available at http://space.mit.edu/ASC/MARX/
>
> Witt, H. J., Mao, S., & Schechter, P. L. 1995, ApJ, 443, 18
>
> Wyithe, J. S. B., Webster, R. L., Turner, E. L., & Mortlock, D. J. 2000, MNRAS, 315, 62
>
> Yonehara, A. 2001, ApJ, 548, L127

| TABLE 1 | | | | |
| --- | --- | --- | --- | --- |
| RESULTS FROM FITS TO THE SPECTRA OF H 1413+117 | | | | |
| Fit | Source1 | Model | Parameter2 | Value3 |
| 1 | H 1413+1174 | PL, neutral absorption | $`\mathrm{\Gamma }_A`$ | 1.8$`{}_{}{}^{+0.8}{}_{0.8}{}^{}`$ |
|  | (Image A) | at source, and a Gaussian | $`N_{\mathrm{H},\mathrm{A}}`$ | (2.3$`\pm 2.2`$) $`\times `$ 1023 cm-2 |
|  |  | emission line at source. | Eline | 6.3$`{}_{}{}^{+0.1}{}_{0.1}{}^{}`$ keV |
|  |  | The fit is performed to the | $`\sigma _{line}`$ | $`<`$ 0.14 keV |
|  |  | spectrum of image A of H1413. | EWline | 1.1$`{}_{}{}^{+0.8}{}_{0.5}{}^{}`$ keV |
|  |  |  | C-statistic/nbins5 | 4.5/13 |
|  |  |  | $`\chi ^2/\nu `$ | 4.0/7 |
|  |  |  | $`P\left(\chi ^2/\nu \right)^6`$ | 0.78 |
| 2 | H 1413+117 | PL, neutral absorption | $`\mathrm{\Gamma }_A`$ | 2.7$`{}_{}{}^{+0.7}{}_{0.7}{}^{}`$ |
|  | (Image A) | at source. | $`N_{\mathrm{H},\mathrm{A}}`$ | (5$`\pm 2`$) $`\times `$ 1023 cm-2 |
|  |  | The fit is performed to | C-statistic/nbins | 19.2/13 |
|  |  | the spectrum of image A of H1413. | $`\chi ^2/\nu `$ | 13.6/10 |
|  |  |  | $`P\left(\chi ^2/\nu \right)`$ | 0.19 |
| 3 | H 1413+117 | PL, neutral absorption | $`\mathrm{\Gamma }_{BCD}`$ | 1.7 f7 |
|  | (Images B + C + D ) | at source, and a narrow Gaussian | $`N_{\mathrm{H},\mathrm{BCD}}`$ | (3$`{}_{}{}^{+3}{}_{2}{}^{}`$) $`\times `$ 1023 cm-2 |
|  |  | emission line at source. | Eline | 6.4$`{}_{}{}^{+0.1}{}_{0.1}{}^{}`$ keV |
|  |  | The fit is performed to the | $`\sigma _{line}`$ | 0.01 keV f |
|  |  | the combined spectrum | EWline | 0.35$`{}_{}{}^{+0.50}{}_{0.35}{}^{}`$ keV |
|  |  | of images B, C, and D. | C-statistic/nbins | 10.64/12 |
|  |  |  | $`\chi ^2/\nu `$ | 7.42/8 |
|  |  |  | $`P\left(\chi ^2/\nu \right)`$ | 0.49 |
| 4 | H 1413+117 | PL, neutral absorption | $`\mathrm{\Gamma }_{BCD}`$ | 1.7$`{}_{}{}^{+0.7}{}_{0.7}{}^{}`$ |
|  | (Images B + C + D) | at source. | $`N_{\mathrm{H},\mathrm{BCD}}`$ | (4$`{}_{}{}^{+3}{}_{2}{}^{}`$) $`\times `$ 1023 cm-2 |
|  |  | The fit is performed to | C-statistic/nbins | 10.64/12 |
|  |  | the combined spectrum | $`\chi ^2/\nu `$ | 9.87/9 |
|  |  | of images B, C, and D. | $`P\left(\chi ^2/\nu \right)`$ | 0.36 |
| 5 | H 1413+117 | PL, neutral absorption | $`\mathrm{\Gamma }`$ | 1.6$`{}_{}{}^{+0.8}{}_{0.9}{}^{}`$ |
|  | Simultaneous fit | at source and a narrow Gaussian | $`N_\mathrm{H}`$ | (2.4$`{}_{}{}^{+2.8}{}_{2.4}{}^{}`$) $`\times `$ 1023 cm-2 |
|  | to images A and | emission line at source. | Eline | 6.3$`{}_{}{}^{+0.1}{}_{0.1}{}^{}`$ keV |
|  | B+C+D | The fit is performed | $`\sigma _{line}`$ | 0.06$`{}_{}{}^{+0.3}{}_{0.06}{}^{}`$ keV |
|  |  | simultaneously to the | EWA | 1.27$`{}_{}{}^{+0.59}{}_{0.59}{}^{}`$ keV |
|  |  | the combined spectrum | EWBCD | 0.20$`{}_{}{}^{+0.36}{}_{0.20}{}^{}`$ keV |
|  |  | of image A and | C-statistic/nbins | 22.5/27 |
|  |  | the combined spectrum | $`\chi ^2/\nu `$ | 18.0/27 |
|  |  | of images B, C, and D. | $`P\left(\chi ^2/\nu \right)`$ | 0.52 |

1All model fits include Galactic absorption towards the source (Stark et al. 1992).
2All absorption-line parameters are calculated in the quasar rest frame.
3All errors are for 90% confidence with all parameters taken to be of interest except absolute normalization. All best fit parameters and associated errors are taken from fits using the C-statistic.
4For results on the combined spectrum of all images of H 1413+117 see also Oshima et al. (2000) and Gallagher et al. (2002).
5 nbins is the number of pulse height analysis (PHA) bins in the spectrum.
6 $`P(\chi ^2/\nu )`$ is the probability of exceeding $`\chi ^2`$ for $`\nu `$ degrees of freedom.
7 The symbol f indicates that the parameter is frozen in the spectral fit.


format: text

text:
# Weighted evolving networks: coupling topology and weights dynamics

## Abstract

We propose a model for the growth of weighted networks that couples the establishment of new edges and vertices and the weights’ dynamical evolution. The model is based on a simple weight-driven dynamics and generates networks exhibiting the statistical properties observed in several real-world systems. In particular, the model yields a non-trivial time evolution of vertices’ properties and scale-free behavior for the weight, strength and degree distributions.

The last years have witnessed a hectic activity devoted to the characterization and understanding of networked structures as diverse as ecological and biological systems or the Internet and the WWW Barabasi:2000 ; Amaral:2000 ; mdbook ; psvbook . These networks generally exhibit complex topological properties such as the small-world phenomenon watts98  and scale-free behavior Barabasi:2000 . The need for explaining these complex topological features has led to a new modeling framework that focuses on the dynamical evolution and growth of networks Barabasi:1999 . In this perspective, a wide array of models aimed at capturing various properties of real networks have been formulated mdbook . These models, however, do generally consider only the topological structure and do not take into account the interaction strength –the weight of the link– that characterizes real networks. For example, the diversity of the predator-prey interaction and of metabolic reactions is considered a critical ingredient of ecosystems Pimm ; Krause:2003  and metabolic networks Almaas:2004 , respectively; in social systems, the weight of interactions is very important in the characterization of the corresponding networks Granovetter . Similarly, the Internet traffic psvbook  or the number of passengers in the airline network Amaral:2000 ; Guimera:2003 ; Barrat:2003  are crucial quantities in the study of these systems. Interestingly, recent studies Almaas:2004 ; Barrat:2003 ; Li:2003a ; Garla:2003  have shown that weighted networks exhibit additional complex properties such as broad distributions and non-trivial correlations of weights that do not find an explanation just in terms of the underlying topological structure.

In this Letter, we define a simple model for weighted networks that considers the basic evolution of the system as driven by the weight properties of edges and vertices. In addition, differently from early models proposed in the past Yook:2001 ; zheng03 , we allow the dynamical evolution of weights during the growth of the system. This mimics the evolution and reinforcements of interactions in natural and infrastructure networks. The generated networks display power-law behavior for the weight, degree and strength distributions with non-trivial exponents depending on the unique parameter defining the model’s dynamics. These results suggest that the inclusion of weights in networks modeling naturally explains the diversity of scale-free behavior empirically observed in real networked structures. Strikingly, the weight-driven growth recovers an effective preferential attachment for the topological properties, providing a microscopic explanation for the ubiquitous presence of this mechanism.

Weighted networks are usually described by a matrix $`w_{ij}`$ specifying the weight on the edge connecting the vertices $`i`$ and $`j`$, with $`i,j=1,\mathrm{},N`$ where $`N`$ is the size of the network ($`w_{ij}=0`$ if the nodes $`i`$ and $`j`$ are not connected). In the following we will consider only the case of symmetric weights $`w_{ij}=w_{ji}`$. Prototypical examples of weighted networks can be found in the world-wide airport network (WAN) Guimera:2003 ; Barrat:2003  and the scientific collaboration network (SCN) newmancoll ; vicsek . In the airport network each given weight $`w_{ij}`$ is the number of available seats on direct flights connections between the airports $`i`$ and $`j`$. In the SCN the nodes are identified with authors and the weight depends on the number of co-authored papers newmancoll ; Barrat:2003 . A first topological characterization of networks is obtained with the vertices degree $`k`$, number of connected neighbors, and the associated probability distribution $`P(k)`$. Complex networks often exhibits a power-law degree distribution $`P(k)k^\gamma `$ with $`2\gamma 3`$ and measurements of weighted networks confirm this evidence Guimera:2003 ; Barrat:2003 . Similarly, a first characterization of weights is obtained by the distribution $`P(w)`$ that any given edge has weight $`w`$. In complex networks such as the WAN and SCN this distribution is heavy tailed and spans several orders of magnitude Barrat:2003 ; Li:2003a . Along with the degree of a node, a very significative measure of the network properties in terms of the actual weights is obtained by looking at the vertex strength $`s_i`$ defined as Yook:2001 ; Barrat:2003

$$s_i=\underset{j𝒱(i)}{}w_{ij},$$
(1)
where the sum runs over the set $`𝒱(i)`$ of neighbors of $`i`$. The strength of a node integrates the information about its connectivity and the weights of its links, and can be considered as the natural generalization of the connectivity. For instance, in the case of the WAN the strength provides the actual traffic going through the vertex $`i`$ and is an obvious measure of the size and importance of each airport. This quantity obviously increases with the vertex degree $`k_i`$ and usually displays the power-law behavior $`sk^\beta `$, with the exponent $`\beta `$ varying as a function of the specific network. Finally, we note that in most cases, the empirical distribution of $`s`$ is heavy-tailed Barrat:2003 , analogously to the power-law decay of the degree distribution.

Previous approaches to the modeling of weighted networks focused on growing topologies where weights were assigned statically, i.e. once for ever, with different rules related to the underlying topology Yook:2001 ; zheng03 . These mechanisms, however, overlook the dynamical evolution of weights according to the topological variations. We can illustrate this point in the case of the airline network. If a new airline connection is created between two airports it will generally provoke a modification of the existing traffic of both airports. In general, it will increase the traffic activity depending on the specific nature of the network and on the local dynamics. In the following we will implement a model that takes into account the coupled evolution in time of topology and weights.

The model dynamics starts from an initial seed of $`N_0`$ vertices connected by links with assigned weight $`w_0`$. At each time step, a new vertex $`n`$ is added with $`m`$ edges that are randomly attached to a previously existing vertex $`i`$ according to the probability distribution

$$\mathrm{\Pi }_{ni}=\frac{s_i}{_js_j}.$$
(2)
This rule relaxes the usual degree preferential attachment, focusing on a strength driven attachment in which new vertices connect more likely to vertices handling larger weights and which are more central in terms of the strength of interactions. The weight of each new edge is fixed to a value $`w_0`$. Moreover, the presence of the new edge $`(n,i)`$ will introduce variations of the existing weights across the network. In particular, we consider the local rearrangements of weights between $`i`$ and its neighbors $`j𝒱(i)`$ according to the simple rule

$$w_{ij}w_{ij}+\mathrm{\Delta }w_{ij},$$
(3)
where

$$\mathrm{\Delta }w_{ij}=\delta \frac{w_{ij}}{s_i}.$$
(4)
This rule considers that the establishment of a new edge of weight $`w_0`$ with the vertex $`i`$ induces a total increase of traffic $`\delta `$ that is proportionally distributed among the edges departing from the vertex according to their weights (see Fig. 1), yielding $`s_is_i+\delta +w_0`$. At this stage, it is worth remarking that while we will focus on the simplest model with $`\delta =const`$, different choices of $`\mathrm{\Delta }w_{ij}`$ with heterogeneous $`\delta _i`$ or depending on the specific properties of each vertex ($`w_{ij},k_i,s_i`$) can be considered Note . Finally, after the weights have been updated the growth process is iterated by introducing a new vertex with the corresponding re-arrangement of weights.

The weight driven attachment (Eq. (2)) appears to be a plausible mechanism in many technological networks. In the Internet new routers connect to more central routers in terms of bandwidth and traffic handling capabilities and in the airport networks new connections are generally established to airports with a large passenger traffic. Even in the SCN this mechanism might play a role since an author with more co-authored papers is more visible and open to further collaborations. Also the weights dynamics can be traced back to simple dynamical processes. Indeed, the quantity $`w_0`$ sets the scale of the weights and without loss of generality we can use the rescaled quantities $`w_{ij}/w_0`$, $`s_i/w_0`$ and $`\delta /w_0`$. For the sake of simplicity we can set $`w_0=1`$ and the model depends only on the dimensionless parameter $`\delta `$, that is the fraction of weight which is ‘induced’ by the new edge onto the others. According to the value of $`\delta `$, different scenarios are possible. If the induced weight is $`\delta w_0=1`$ we mimic situations in which an appreciable fraction of traffic generated by the new connection will be dispatched in the already existing connections. This is plausible in the airport networks where the transit traffic is rather relevant in hubs. In the case of $`\delta <1`$ we face situations such as the SCN where it is reasonable to consider that the birth of a new collaboration (co-authorship) is not triggering a more intense activity on previous collaborations. Finally, $`\delta >1`$ is an extreme case in which a new edge generates a sort of multiplicative effect that is bursting the weight or traffic on neighbors.

In the model, the time is measured with respect to the number of vertices added to the graph, resulting in the definition $`t=NN_0`$, and the natural time scale of the model dynamics is the network size $`N`$. The network’s evolution can be inspected analytically by studying the time evolution of the average value of $`s_i(t)`$ and $`k_i(t)`$ of the $`i`$-th vertex at time $`t`$, and by relying on the continuous approximation that treats $`k`$, $`s`$ and the time $`t`$ as continuous variables Barabasi:2000 ; mdbook . When a new edge is added to the network, the strength $`s_i`$ of vertex $`i`$ can increase either if the edge connects directly to $`i`$ or to one of its neighbor. The evolution equation for $`s_i`$ is thus given by

$`{\displaystyle \frac{ds_i}{dt}}`$ $`=`$ $`m{\displaystyle \frac{s_i}{_js_j}}(1+\delta )+{\displaystyle \underset{j𝒱(i)}{}}m{\displaystyle \frac{s_j}{_ls_l}}\delta {\displaystyle \frac{w_{ij}}{s_j}}`$ (5)
$`=`$ $`{\displaystyle \frac{2\delta +1}{2\delta +2}}{\displaystyle \frac{s_i(t)}{t}},`$
where the last expression is recovered by noticing that each added link increases the total strength by an amount equal to $`2+2\delta `$ thus implying that $`_is_i(t)2m(1+\delta )t`$.

This equation can be readily integrated with the initial condition $`s_i(t=i)=m`$, yielding

$$s_i(t)=m\left(\frac{t}{i}\right)^{\frac{2\delta +1}{2\delta +2}}$$
(6)
Similarly, the rate equation for the degree evolution reads

$$\frac{dk_i}{dt}=m\frac{s_i(t)}{_js_j(t)}=\frac{s_i(t)}{2(1+\delta )t},$$
(7)
that implies $`k_i(t)=s_i(t)/(2\delta +1)`$. The proportionality relation $`sk`$, that implies $`\beta =1`$, is particularly relevant since it states that the weight-driven dynamics generates in Eq. (2) an effective degree preferential attachment that is parameter independent. This highlights an alternative microscopic mechanism accounting for the presence of the preferential attachment dynamics in growing networks.

The behavior of the various statistical distributions can be computed by noticing that the time at which the node $`i`$ enters the network is uniformly distributed in \[0,t\] and

$$P(k,t)=\frac{1}{t+N_0}_0^t\delta (kk_i(t))𝑑i,$$
(8)
where $`\delta (x)`$ is the Dirac delta function. Solving the above integral and considering the infinite size limit $`t\mathrm{}`$ yields $`P(k)k^\gamma `$ with

$$\gamma =\frac{4\delta +3}{2\delta +1}.$$
(9)
It is straightforward to obtain the same behavior $`P(s)s^\gamma `$ for the strength distribution since $`s`$ and $`k`$ are proportional. This result shows that the obtained graph is a scale-free network described by an exponent $`\gamma [2,3]`$ that depends on the value of the parameter $`\delta `$. In particular, when the addition of a new edge doesn’t affect the existing weights ($`\delta =0`$), the model is topologically equivalent to the Barabasi-Albert model Barabasi:1999  and the value $`\gamma =3`$ is recovered. For larger values of $`\delta `$ the distribution is progressively broader with $`\gamma 2`$ when $`\delta \mathrm{}`$. This indicates that the weight-driven growth generates scale-free networks with exponents varying in the range of values usually observed in the empirical analysis of networked structures Barabasi:2000 ; mdbook ; psvbook . Noticeably the exponents are non-universal and depend only on the parameter $`\delta `$ governing the microscopic dynamics of weights. The model therefore proposes a general mechanism for the occurrence of varying power-law behaviors without resorting on more complicate topological rules and variations of the basic preferential attachment mechanism.

In order to check the analytical predictions we performed numerical simulations of networks generated by using the present model with different value of $`\delta `$, minimum degree $`m`$ and varying size $`N`$. In Fig. 2 we show the behavior of the vertices’ strength versus time for different value of $`\delta `$, recovering the behavior predicted analytically. We also report the average strength $`s_i`$ of vertices with degree $`k_i`$ and confirm that $`\beta =1`$ as well as the value of the pre-factor $`2\delta +1`$. Finally, in Fig. 3A we plot the probability distribution $`P(s)`$ that exhibits a power-law behavior in excellent agreement with the theoretical predictions.

Similarly to the previous quantities, it is possible to obtain analytical expressions for the evolution of weights and the relative statistical distribution. The weight $`w_{ij}`$ increases by the addition of a new link either on $`i`$ or on $`j`$ and the corresponding rate equation can be written as

$$\frac{dw_{ij}}{dt}=m\frac{s_i}{_js_j}\delta \frac{w_{ij}}{s_i}+m\frac{s_j}{_js_j}\delta \frac{w_{ij}}{s_j}.$$
(10)
This equation can be solved with the initial condition $`w_{ij}(t_{ij})=1`$, where $`t_{ij}=`$max$`(i,j)`$ is the time at which the edge is established, yielding $`w_{ij}(t)=(t/t_{ij})^{\delta /(\delta +1)}`$ From this expression it is possible to calculate similarly to $`P(k)`$ the probability distribution $`P(w)`$ which is in this case also a power-law $`P(w)w^\alpha `$ where

$$\alpha =2+\frac{1}{\delta }.$$
(11)
The exponent $`\alpha `$ has large variations as a function of the parameter $`\delta `$ and $`P(w)`$ moves from a delta function for $`\delta =0`$ to a very slow decaying power-law with $`\alpha =2`$ if $`\delta \mathrm{}`$. This feature clearly shows that the weight distribution is extremely sensitive to changes in the microscopic dynamics ruling the network’s growth. These predictions can be confirmed by numerical simulations and the results are shown in Figs. 3B and 4.

The model we have presented is possibly the simplest one in the class of weight-driven growing networks. A novel feature in the model is the weight dynamical evolution occurring when new vertices and edges are introduced in the system. This simple mechanism produces a wide variety of complex and scale-free behavior depending on the physical parameter $`\delta `$ that controls the local microscopic dynamics. While a constant parameter $`\delta `$ is enough to produce a wealth of interesting network properties, a natural generalization of the model consists in considering $`\delta `$ as a function of the vertices degree or strength. For example, a non-linear generalization of eq. (4) can yield non-linear correlations between strengths and degrees ($`\beta 1`$Note . Similarly, more complicated variations of the microscopic rules may be implemented to mimic in a detailed fashion particular networked systems. In this perspective the present model appears as a general starting point for the realistic modeling of complex weighted networks.

###### Acknowledgements.

We thank R. Pastor-Satorras and M. Vergassola for useful discussions. A.B and A. V. are partially funded by the EC-Fet project COSIN IST-2001-33555.


text:
# Emergence of Conduction Channels in Lithium Silicate

## I Introduction

The presence of conduction pathways has been proposed to rationalize the fast conduction in ion conductors like alkali silicates. Based on EXAFS data on the coordination environments in these systems, Greaves proposed the modified random network model for glassesGreaves (1985). It predicts the aggregation of modifier ions into channels lined by non-bridging oxygen atoms. Later NMR studiesYap et al. (1995); Gee and Eckert (1996) reported both inhomogeneous and homogeneous distributions of cations in the structure with a clustering of alkali ions in the low concentration regime. These results are backed by molecular dynamics (MD) simulations of alkali silicate where also only for low alkali concentrations a mild clustering of alkali ions has been found Heuer et al. (2002); Horbach et al. (2001).

In his analysis of MD simulations of sodium silicate, Jund et al Jund et al. (2002) focussed on the most mobile cations. They introduced $`n_{\mathrm{diff}}`$, the number of different ions visiting a given subvolume, to determine those parts of the system taking part in effective long range transport. They report that the regions of the silicate network, visited by the largest number of different ions form a network of blobs and connecting channels (denoted conduction channels). Stated differently, these regions show a significant spatial correlation beyond the nearest-neighbor shell. In particular, these spatial correlations are much stronger than any possible correlation of alkali positions themselves. Thus these conduction channels do not simply emerge from clustering of alkalisJund et al. (2002) .

Interestingly, all static structure factors for sodium-silicate exhibit a prepeak at $`q_10.95`$ Å-1, corresponding to the next-nearest Na-Na or Si-Na neighbors, as shown experimentallyMeyer et al. (2002) and numericallyHorbach et al. (2002); Sunyer et al. (2003a). Thus the structure displays some long-range correlations. One may be tempted to relate these long-range correlations to the channels discussed above. The validity of this conclusion, however, is not obvious, since the structure refers to all ions whereas the channels emerge from the properties of the specific subset of regions, visited by many different ions.

In this paper we analyze the question to which degree these conduction channels are predetermined by the structure. This question is central because it elucidates the role of the network for long-range transport. Two extreme scenarios are conceivable. First, following the concepts of Greaves, the network may supply channels where ion conduction is strongly favored due to the aggregation of the relevant structural elements (in Greaves work: non-bridging oxygens) which favor fast dynamics. In this scenario these relevant structural elements would also display long-range spatial correlations. Second, the structure might only determine the local mobility of an ion such that the relevant structural elements do not display these long-range correlations. Rather the formation of conduction channels is a simple statistical process in analogy to the formation of percolation paths in random-energy or random-barrier lattice models of ion conductionDyre and Schrøder (2000).

In previous work we have presented a method how to identify individual lithium sites which are basically time-invariant below the glass transitionLammert et al. (2003). In particular, we could show that the continuous trajectories of the lithium ions can be interpreted as hops between these sites. With this information at hand we can identify the conduction channels by counting the number of different ions visiting a specific site during a long MD simulation. This discretization of the system is somewhat more adapted to the conduction process than the simple tiling by small cubes as done by Jund et al Jund et al. (2002). In particular, we are very sensitive to the effect of structural elements like bridging or non-bridging oxygens on the alkali dynamics.

The paper is organized as follows. In Section II we discuss the technical aspects of our work. Section III contains the results which are finally discussed in Section IV.

## II Technical aspects

The data for this investigation are based on MD simulations of lithium silicate $`(Li_2O)(SiO_2)`$. The system containing 1152 atoms was propagated under NVT conditions at 640K with 2 fs stepsize for 10 ns simulation time. As a basis for investigations of the cation dynamics we first located the individual sites available to the lithium ions. For this purpose we have discretized the system in very small cubes (0.3 Å) 3 and determined the occupation number in every cube. This is a simple way to sample the effective potential energy landscape provided by the silicate network. In order to distinguish regions which indeed belong to ionic sites and those which only serve as transition paths between these sites we have used some reasonable cutoff criterion for the identification of those cubes which are part of ionic sites. Subsequent cluster analysis of the remaining cubes finally defines the ionic sites. For our system we have obained 378 different sites. Details about the method and about the properties of the sites are given in an earlier paperLammert et al. (2003).

We characterize the mobility of a site in two different ways. First, in analogy to the work by Jund et al. the number of different ions visiting a site has been determined. It is denoted $`n_{\mathrm{diff}}`$. Sites for which the value of $`n_{\mathrm{diff}}`$ belongs to the 10% highest are denoted A-sites. In contrast, those with the 10% smallest values of $`n_{\mathrm{diff}}`$ are denoted $`\overline{A}`$-sites. Second, we determined the average residence time in a site, $`\tau `$. Analoguosly, sites with the 10% shortest residence times are denoted B-sites, those with the longest 10% residence times $`\overline{B}`$-sites. In general terms, A- and B-sites indicate positions in the network where ions are fast, $`\overline{A}`$\- and $`\overline{B}`$-sites where ions are slow. In this investigation both quantities will be used together as different measures of site mobility.

Furthermore we performed straightforward Monte Carlo simulations for hopping particles on a 2D square lattice model with 1600 sites and periodic boundary conditions. We introduced random site energy disorder and additional barrier disorder among adjacent sites. Both distributions were chosen as constant with energies between 0 and 8 $`k_BT`$. 90% of all sites were populated and the hopping particles interacted via simple excluded volume interaction. The energy disorder was kept constant during the whole run.

## III Results

### III.1 General properties of sites

The two definitions of site mobility are compared in Fig. 1, which shows $`n_{\mathrm{diff}}`$ plotted against $`\tau `$ for all sites. The highest possible value of $`n_{\mathrm{diff}}`$ for a site is basically given by the number of residences permitted by their average duration $`\tau `$: $`n_{\mathrm{diff}}<=t/\tau `$. This implies that a sufficiently short average residence time is a necessary condition for a large number of different ions visiting a site. A broad range of lower values for $`n_{\mathrm{diff}}`$ is also found, caused by multiple visits of identical ions. The correlation between $`\tau `$ and $`n_{\mathrm{diff}}`$ is therefore only weak.

General information on the structural environment of the sites can be obtained from the radial distribution function g(r) of network atoms around lithium. It is shown in Fig. 2 for bridging oxygens (BOs), non-bridging oxygens (NBOs), and silicon atoms. BOs and NBOs have been distinguished according to the number of silicon atoms within 2.4 Å, the distance of the first minimum in the g(r) of silicon and oxygen. One finds that lithium atoms are surrounded by oxygens. Integration up to the first minimum gives average coordination numbers of 4.0 NBOs and 1.2 BOs. In contrast, if the average number of 5.2 oxygen neighbors were divided according to the relative numbers of BOs and NBOs in the system one would obtain 3.3 NBOs and 1.8 BOs. As expected, NBOs are favored as coordination partners of lithium. Greaves (1985); Kamijo et al. (1996); Florian et al. (1996); Sunyer et al. (2003b).

### III.2 Spatial correlations of sites

Now we analyze the question to which degree the mobilities of the different sites are correlated. The presence of channels, defined as regions of the system where many different ions pass would imply that A-sites (and correspondingly $`\overline{A}`$-sites) are spatially correlated over larger distances. Furthermore, also possible spatial correlations of B-sites are analyzed.

It can be expected that correlations on the length scale of the nearest-neighbor distance occur due to trivial reasons. For a fast B-site it can be expected that there is a low barrier to an adjacent site. Accordingly, also the residence time in this adjacent site can be expected to be small because a fast backjump to the original site would be possible. This trivial correlation is limited to next neighbors. If B-sites were arranged into channels, they must show an additional correlation over longer distances. The presence of a continuous pathway of B-sites through the system would result in an increased probability to find other B-sites among the second neighbors of a given B-site.

Apart from A- and $`\overline{A}`$-sites (and analogously B- and $`\overline{B}`$-sites) we also selected a subset of 10% of sites with intermediate $`\tau `$ and $`n_{\mathrm{diff}}`$, respectively. Each of the three groups contained 37 sites, i.e. ca. 10% of the total number. Because of the low number of sites we did not compute their complete pair distribution functions g(r). Instead we counted directly the number of neighbors from each of the three groups around a given central site. The limits of the first and second neighbor shell were taken from the minima in the g(r) of all sites at 4.2 Å and 6.8 Å. The results for slow, intermediate and fast central sites are compared in Fig. 3 and Fig. 4 for the first and for the second shell respectively.

All values are given as a fraction of the total number of neighbors. The statistical value of $`37/378=9.8\%`$ for the contribution of each group is shown by the dashed lines. Generally, there is a clear tendency that a site will be surrounded by first neighbors of similar mobility. The strongest correlation is found among the fast sites. While the findings are similar for both $`\tau `$ and $`n_{\mathrm{diff}}`$, they are more pronounced in the latter case. For the second shell of neighbors, the correlations almost vanish for the residence time $`\tau `$, i.e. B-sites (and $`\overline{B}`$-sites) are spatially not correlated. In contrast, as already reported by Jund et al. A-sites are significantly correlated beyond the nearest-neighbor shell. The correlation for $`\overline{A}`$-sites is much weaker, but still present.

### III.3 Spatial correlations of relevant network properties

Next we turn to the question in which sense the conduction channels are predetermined by the structure. To answer this question we have to proceed in two steps. First, we have to check which property of the network significantly favors highly mobile sites. As will be shown below it is the presence of bridging oxygens (BOs) which displays the strongest correlation with the mobility of a nearby site. Second, possible long-range spatial correlations of this network property (here: occurrence of BOs) have to be analyzed. More specifically we check whether sites with a high/low number of BOs are spatially correlated. Again we have to be aware that trivial nearest-neighbor correlations will be present. If some site displays an increased number of BOs as compared to the average site it can be expected that a small number of neighbor sites shows the same correlation because one BO typically belongs to 2 lithium ions.

Using the sets of sites already defined above, we first count the number of BOs, NBOs and Si atoms up to the 1st minimum distances of 2.9Å, 3.0Å, and 4.1 Å respectively. The results are shown in Fig. 5. As expected from the above-mentioned analysis of coordination numbers, sites are on average surrounded by $``$ 4 NBOs, while only 0.5 to 1.5 BOs are found. For comparison, lines are shown for 1.8 BOs and 3.3 NBOs, i.e. the statistical values (see above). Interestingly, the strongest influence on the mobility of a site is shown for the number of coordinated BOs, which decreases by $`1`$ from slow to fast sites. The number of silicon neighbors shows only a weak decrease in the same direction. For NBOs the correlation is reverse, but much weaker even on an absolute scale. The results are very similar when using $`\tau `$ or $`n_{\mathrm{diff}}`$ to determine the mobility of the sites.

This result implies that in most cases a single BO in the nearest-neighbor shell is sufficient to immobilize ions on the sites around the BO. A possible underlying structure in the network that could determine the layout of conduction channels should therefore be related to the distribution of BOs. Similarly to the method described above, we sorted all sites into three groups according to the number of BOs in their 1st shell closer than 2.9 Å. Then the fraction of neighbors from each of these groups was counted. The results are given in Fig. 6. The dashed line again marks the statistical value, which is 0.33 in this case, corresponding to groups of one third of all sites. The expected trivial correlation among nearest neighbors is, of course, found. Most important, no correlation to the second nearest neighbor shell is found. We checked that also the NBOs do not display any correlations beyond the nearest neighbor shell. The lack of spatial correlations beyond the nearest neighbors strongly suggests that even the distribution of bridging oxygen atoms as the most relevant network property does not determine complete conduction channels in the system.

## IV Discussion

The presence of only one BO drastically slows down the ionic transport through nearby sites. An indication how this slowdown is effected can be obtained from the two examples shown in Fig.7. Here the local network structure is shown around one slow site (left) and one fast site. Cations are shown at the center as solid balls. The surrounding silicate network is shown as tubes, . The site as determined by usLammert et al. (2003) is represented by the translucent cloud surrounding the cation. Connections are drawn between the lithium ion and all oxygen atoms within the first neighbor shell. For the fast site, all six coordinating oxygen atoms are NBOs, which are bound to four different silicon atoms. For the slow site there is one BO among the neighbors. This BO is part of a silica chain segment that is also carrying three of the NBOs coordinating the cation. These NBOs are bound to three silicon atoms on different sides of the BO, with an additional BO inserted on one side in between. In total, a chain segment of 5 atoms is linked to the cation via four coordinating oxygen atoms. The chain is thus closely wrapped around a large part of the site. This situation suggests that BOs introduce high transport barriers for nearby sites because the silica chains joined by them block a large fraction of possible jump pathways.

The pathways through the system used for effective ionic transport follow continuous channels of sites visited by many ions. The passage of a large number of different ions through a specific site depends on two conditions. (i) A comparatively low average residence time of the sites is necessary but not sufficient to enable the passage of many ions. (ii) The presence of coordinating BOs strongly decreases the number of ions reaching a site, probably by introducing a local barrier. Both the distribution of sites with low $`\tau `$ and that of BOs in the structure show no long range order comparable to that found in the pathways traced by sites with high $`n_{\mathrm{diff}}`$. They can therefore only impose constraints on the layout of these pathways, but not determine them. Therefore the emergence of conduction channels does not seem to require preformed channels provided by the network structure.

Now we show that on a qualitative level similar effects can be observed for the simple lattice model, introduced in Section II. After performing long Monte Carlo simulations we determined the A- and B-sites similarly to our procedure for the silicate system. They are marked by filled dots and open circles, respectively. Interestingly, the three main observations from the silicate system are reproduced: (i) By definition the structure (here: the local energetic parameters) does not display correlations beyond the nearest neighbor shell. (ii) The B-sites are not spatially correlated beyond the nearest-neighbor shell. Actually, one clearly sees that sites with short residence time typically occur as pairs, because of the trivial correlations discussed above. (iii) The A-sites are spatially correlated, i.e. they are concentrated in a few larger blobs. These regions contain some of the groups of sites with low $`\tau `$, but also some sites with higher residence times. Correspondingly, some sites with low $`\tau `$ are not incorporated into regions with high $`n_{\mathrm{diff}}`$. Sites visited by many different ions can of course not be completely isolated. The ions have to move away through other sites in the neighborhood to allow the following ions to enter. In our lattice model the regions best suited for this process are located by chance. The ions benefit from sites with very low residence times, they must reach them via other sites that are also not too slow, and they must avoid high barriers distributed through the lattice. Interestingly, the pathways are somewhat extended in both dimensions. This suggests that regions where many different particles may diffuse are extended in all directions. Actually, the channels, as obtained from the molecular dynamics simulationsJund et al. (2002) show this property as well.

In summary, we have analyzed the nature of the conduction pathways in lithium silicate as well as a simple random lattice model. In both cases we do find spatial correlations in regions of high mobility. These structures do not seem to be supplemented by corresponding spatial correlations of network properties, following from the lack of long-range correlations of sites with BOs. Thus we would like to suggest that the formation of conduction channels is rather a statistical selection process in analogy to the formation of percolation paths in simple lattice models.

###### Acknowledgements.

We thank M. Korth and H. Krieg for the help with the MC-simulations. Furthermore we acknowledge very helpful conversations with P. Jund. H. L. also acknowledges the support by the Fonds der Chemischen Industrie.


text:
# Feasibility of a magnetic suspension for second generation Gravitational Wave interferometers

## 1 Introduction

Magnetic levitation systems find several applications  in physics and engineering, ranging from small force measurements to transportation systems. In this paper, we analyze the possibility of using magnetic levitation as a future upgrade of the sophisticated suspensions used in long baseline interferometric Gravitational Wave (GW) detectors like, for example, Virgo .

The detection of GW is one of the most challenging fields of contemporary physics  and interferometers appear to be the most promising detectors looking for the first direct evidence of GW. Presently, several long baseline antennas have already started operation , or are still under construction . Due to the extreme weakness of the GW signals, the main requirement for these detectors is the isolation of the optical components – used as free-falling test masses – from any kind of external disturbance. For this reason, the whole optical path (e.g. three km per arm for Virgo) is under ultra high vacuum and all the mirrors are suspended to complex vibration isolation systems, which strongly reduce the effect of seismic vibration, that is the main limiting noise below few Hz.

The idea of using stable levitation for the mirror suspensions of long baseline GW detectors has been proposed by Drever  for magnetic levitation, and Giazotto  for electrostatic levitation. The main advantage of these systems is the possibility to hold the mirror without direct contact; therefore, these configurations avoid some noise contributions due to the presence of suspension wires – for example the thermal noise associated with violin modes.

In the case of magnetic levitation, the presence of big magnets directly attached to the mirrors would decrease the mechanical quality factor of the internal modes, thus increasing the associated thermal noise, and should make the mirror position directly sensitive to external electro-magnetic (EM) noise. On the other hand, using an electrostatic levitation system is safe from the point of view of the mirror quality factor.

To obtain a stable levitation, it is necessary to use a servo-loop to control the vertical position of the mirror. In this case, the noise of the position sensor should be low enough, in order not to spoil the sensitivity of the antenna. For example, for the Virgo interferometer a sensitivity of $`10^{17}\mathrm{m}/\sqrt{\mathrm{Hz}}`$ at 10 Hz is expected. Assuming a (optimistic) vertical to horizontal coupling factor of 1%, the maximum acceptable noise for the sensor is $`10^{15}\mathrm{m}/\sqrt{\mathrm{Hz}}`$. This specification is quite severe, although not impossible, and requires an interferometric readout system.

All these questions must be carefully investigated before concluding that direct levitation of the mirrors can be actually applied to GW detectors. Another possibility, studied in this paper, is the usage of a levitation system as an intermediate stage of a Virgo-like multi-stage seismic attenuator.

The Virgo seismic suspension, called Superattenuator (SA), is a pendulum with blade springs at each of its five stages, which provides a vibration attenuation in six degrees of freedom by more than 200 dB in the whole measurement frequency band of the antenna (4 Hz $`÷`$ 10 kHz). A detailed description of the SA can be found in ; in the following, we will only review its most relevant features. Fig. 1 shows a Virgo SA which can be conceptually divided into two parts.

* The upper one, including the 6-meter inverted pendulum supporting the five attenuation stages.
* The lower one, called payload, made of an intermediate mass – the marionetta –, supporting the mirror itself, and of the reference mass, a cylindrical shell suspended around the mirror with same mass and same barycenter.

Tab. 1 lists the masses of the different SA components. The mirror orientation and position can be controlled using coils as actuators. They are either attached to the last stage of the upper part of the SA and acting on magnets glued on the marionetta, or fixed to the reference mass and acting on small magnets directly fixed on the mirror edges. The lower part of the SA is in ultra-high vacuum ($`<10^7`$ mbar) to avoid contamination of the mirror surfaces by hydrocarbonates in particular, while the upper part is only under high vacuum ($`10^6`$ mbar), because of the outgassing of several components of the suspension. To avoid the contamination of the mirrors, the two vacuum chambers are separated by a roof with a small conductance hole for the suspension wire.

Due to the low frequency thermal drifts of the system, the centering of the wire in conductance hole is quite difficult and requires periodic re-adjustment, or the use of a closed loop control for long term operation.

A possible alternative solution of the problem is the levitation of the marionetta, as shown in Fig. 2. In this way, the upper and lower parts of the vacuum chamber could be completely isolated by a dielectric and transparent roof, eliminating the conductance hole. This configuration is much less sensitive to the noise problems described earlier, as both the sensor noise and the external EM noise are injected at the level of the marionetta and are thus filtered by the last stage of the suspension, providing an additional attenuation of more than two orders of magnitude at 10 Hz. Indeed, the specification on the sensor noise can be relaxed to $`10^{13}\mathrm{m}/\sqrt{\mathrm{Hz}}`$.

The main goal of this work is thus to evaluate the possibility of suspending with magnetic levitation a given mass (on the order of the sum of the masses of the mirror and of the reference mass) at manageable cost.

In the following sections, we analyze the feasibility of such a levitated system. In Section 2, the general principle of stable levitation  is shortly recalled. Then, the numerical procedures used to compute the forces between the different magnets are described. In Section 2.2, the simulation results for magnetic levitation applied to the Virgo suspensions are presented. An experimental test – described in Section 2.3 – shows the feasibility of the magnetic levitation; the experimental results obtained on a small levitated suspension are compared with the predictions of some numerical simulations in order to check the validity of the model.

## 2 Magnetic suspension system

As it is well known, a stable levitation in a constant magnetic field can be obtained only with superconducting and diamagnetic materials . Using superconducting materials requires low temperature operation; thus, their integration in the seismic suspension of GW detectors looks complex and expensive. Moreover, for diamagnetic materials, levitating very small pieces requires intense fields, which makes the suspension of a payload of several tens of kilograms (see Tab. 1) almost impossible. On the other hand, using a magnetic field of variable intensity controlled by a feedback system in association with permanent magnets and electromagnets is more promising.

The general principle of such device can be seen in Fig. 2 and 3: a permanent magnet of mass $`m_2`$ is attached to the mass to be levitated (for the Virgo SA, the payload and the mirror) and is attracted by a fixed magnet of mass $`m_1`$ (on the SA upper part) which exactly balances the gravitational force. In this configuration, the equilibrium position is stable in the horizontal direction, but unstable in the vertical one. This configuration can become stable if a coil acting on the levitated magnet is added, with a current intensity $`I`$ depending on the vertical position of the levitated mass. This position can be measured with any type of position sensor (a shadow meter for our experiment, see Section 2.3). In principle, the fixed magnet $`m_1`$ could be removed, but in this case, a large DC current would be necessary to balance the weight of the levitated piece.

### 2.1 The numerical simulation

Testing the feasibility of such technique is not enough, as one cannot simply scale a system made of point-like magnets to a large device using extended magnets, such as the one needed for a Virgo SA: the dipole approximation is no more valid to compute the magnetic force. Therefore, an accurate calculation of the forces between the fixed and the levitated magnets is needed.

Numerically, the magnetic field generated by a couple of big magnets is computed by dividing the large pieces into infinitesimal volumes, as shown in Fig. 4. Using the corresponding magnetization per unit volume, it is possible to apply the dipole approximation between any two such small regions (one in each magnet), provided that their separation is much larger than their sizes. The final force is the sum of all these infinitesimal contributions. For instance, in Fig.4, the vector $`\stackrel{}{r}`$ gives the position of an infinitesimal volume of the fixed magnet (with respect to the origin $`O`$ of the reference system, chosen at the center of the fixed magnet), while $`\stackrel{}{r^{}}`$ points on the levitated magnet. $`\mathrm{\Delta }\stackrel{}{r}`$ is defined as the separation vector between the two considered volumes.

The computation of the forces between extended magnets has been implemented in MATLAB . The main aim of this tool is to study the best configuration, by changing the dimensions of the two magnets and their separation. Of course, the accuracy of the computation depends on the dimensions of the infinitesimal volume adopted for the simulation. In each computation, they have been reduced iteratively until convergence (within 2%). In our numerical code, the shape, the size and the separation of the magnets are free parameters; for simplicity, we use parallelepipedal magnets.

### 2.2 Simulation results

For a given geometrical configuration, the force between the magnets depends on the residual magnetization $`B_r`$, a magnet proper parameter. For our computation, we use Nd-Fe-B magnets which have currently the higher residual magnetization ($`B_r=1.3`$ T); the density of this material is $`\rho =7.4\mathrm{g}/\mathrm{cm}^3`$. For each magnet, the dimensions to be optimized are the length, that is the size of the rectangular parallelepiped in the direction the optical axis of the suspended mirror, the width, defined as the size along the transfer dimension and the vertical thickness.

The goal of our study is to determine how the static force between the fixed and the levitated magnet changes with their sizes and their separation. In particular, we use as parameter the ‘free-gap’ $`d`$, i.e. the distance between the two magnets. In this way, we can evaluate the maximum mass that can be levitated in each configuration. In all cases, we need to subtract to the total levitated mass $`M`$ the weight of the levitated magnet $`m`$, so that we get the effective payload (mirror + marionetta) mass $`M_p`$.

Another subject studied here is the dependence of the mass to be suspended and of the horizontal restoring force on the length and the width of the magnets. This allows us to study the performance of the levitated system as a seismic isolator .

The first simulation aims at evaluating the maximum weight of the levitated mass $`M_p`$ as a function of the levitated magnet thickness h. Fig. 5 shows the result: the simulation has been done for two different “free-gap” values ($`d=5`$ cm and $`d=10`$ cm) and for two different upper magnet configurations: $`(20\times 20\times 10)\mathrm{cm}^3`$ and $`(40\times 20\times 10)\mathrm{cm}^3`$. The levitated magnet dimensions are $`(15\times 15\times h)\mathrm{cm}^3`$.

For instance, when the thickness of the levitated magnet to be $`h=2`$ cm, it is possible to suspend about $`50`$ kg when the free-gap is $`10cm`$ and for upper magnet dimensions of $`(20\times 20\times 10)\mathrm{cm}^3`$. As we can see in Fig. 5, the maximum mass to be suspended increases with the thickness and of course strongly depends on the free-gap. An interesting point, apparently counter-intuitive, is that for given values of $`h`$ and $`d`$, the force is smaller for the upper magnet with the larger lateral size. This aspect will be clarified later.

The second step is to evaluate the dependence of $`M_p`$ on the upper magnet thickness H. The result is shown in Fig. 6: simulations have been done for two free-gap values already considered and for two different levitated magnet configurations: the first one, with dimensions $`(15\times 15\times 2)\mathrm{cm}^3`$ and mass $`M_m3.3`$ kg and the second one, with dimensions $`(10\times 10\times 2)\mathrm{cm}^3`$ and mass $`M_m1.5`$ kg; the upper magnet dimensions are $`(20\times 20\times H)\mathrm{cm}^3`$. We can see that it is possible to obtain a levitated mass of $`120.5`$ kg in the configuration $`(15\times 15\times 2)\mathrm{cm}^3`$ and with a gap of 5 cm: as a cross-check, this value corresponds to the one found in Fig. 5 with the same configuration.

Looking at Fig. 6, we can see that it is almost useless to increase the thickness $`H`$ above $`30÷40`$ cm because the force becomes almost constant. As expected, there is still a strong dependence on the gap; for example, we obtain almost the same force with a much smaller suspended magnet, by reducing $`d`$ from 10 cm to 5 cm.

The third simulation studies the evolution of the levitated mass $`M_p`$ when the free-gap $`d`$ varies, as Fig. 7 shows; the upper magnet dimensions are, in this case, $`(20\times 20\times 10)\mathrm{cm}^3`$ and the levitated ones $`(15\times 15\times 2)\mathrm{cm}^3`$. Of course, the suspendable mass value decreases with the increase of the free-gap.

The variation of the suspendable mass $`M_p`$ versus the length L of the upper magnet is the topic of the fourth simulation. As shown in Fig. 8, the calculation has been done for a set of different free-gap values ranging from $`d=5`$ cm to $`d=20`$ cm. The dimensions of the magnets are $`(L\times 20\times 10)\mathrm{cm}^3`$ for the upper one and $`(15\times 15\times 2)\mathrm{cm}^3`$ for the levitated one. As we can see from the Figure, using an upper magnet too long with respect to the free-gap dimension is not an advantage, because the force crosses a maximum and then decreases asymptotically to a constant value. This effect has been already observed in Fig 5. It can be easily explained by noting that the force along the magnetic dipole direction changes its sign when the transverse distance of the two dipoles becomes much larger than their longitudinal one.

The last investigated point is the dependence of the restoring force along the mirror optical axis direction on the variation of the upper magnet length L. This calculation has been done for different free-gap values, $`d=5`$ cm and $`d=20`$ cm, for a levitated magnet $`(15\times 15\times 2)\mathrm{cm}^3`$ and for a constant ‘misalignment’ of 1 mm along the length direction, as shown in Fig. 9. As we can see, the restoring force is always vanishing for $`L\mathrm{}`$; in this case, as proposed by Drever , the system behaves like a pendulum with a very low resonant frequency, i.e. it is in principle a very good seismic isolator.

The interesting point is that, for some gap-length configurations, the restoring force can become negative: in this case the system gets unstable in the horizontal direction. As a consequence, there is a finite length giving a restoring force equal to zero, ranging $`40÷50`$ cm. So, for $`L`$ approaching this value, we can get a stable configuration with very small restoring force even with a finite magnet. For example, if we take a magnet with $`L=40`$ cm, with $`d=5`$ cm we can levitate up to $`95`$ kg (see Fig. 9) with an horizontal restoring force of only 200 N/m, corresponding to a resonant frequency of $`0.23`$ Hz: this is the equivalent of a 45 m-long pendulum. As a comparison, the SA main resonance frequency is around 30 mHz which is equivalent to a 275 m-long simple pendulum.

### 2.3 Experimental test

A magnetic suspension prototype has been realized in Naples  to check the correctness of our models and to verify that we take into account all the relevant effects. Our set-up is sketched in Fig. 3; we use Sm-Co cylindrical magnets which have a residual magnetization $`B_r=0.8`$ T and a density $`\rho =8.3\mathrm{g}/\mathrm{cm}^3`$. The radius is $`R=7`$ mm while the thickness is $`h=12`$ mm for the fixed magnet and 8 mm for the levitated one with a mass $`m`$ of about 10 g. To measure the vertical position of the payload, we use a shadow-meter sensor made of a laser diode partially intersected by the payload lower edge and a photodetector. To get a stable position we feedback on a coil acting on the levitated magnet with a force/current characteristic of about 1 N/A. In this way we are able to hold constant the vertical position of the payload respect to the ground. Using this configuration it is possible to suspend a payload of about 45 g with a free-gap ranging between 1 and 3 cm. The ‘fixed’ magnet is mounted on a micrometric translator allowing to change the free-gap and consequently the vertical force between the magnets. When this force exactly balanced the total weight of the suspended body, the DC current flowing in the coil is zero. The changing of the vertical distance of the magnets with the micrometer results in a non-zero DC current because the force exerted by the coil must balance the difference between the weight and magnets forces. In this way, the current flowing in the coil provides a measurement of the force between the magnets (subtracted of the weight of the suspended mass). Fig. 10 shows the dependence of the measured force between the magnets on the distance between the centers of the magnets. The theoretical curve computed with our simulation model (dashed line) is superimposed to the experimental points. As we can see, the experimental points are in agreement with the model. For comparison, we add the force computed by approximating the extended magnets as point-like dipoles placed in their centers; as expected, the point-like dipole approximation becomes unsatisfactory as the distance decreases and disagrees with the experiment for a distance below 3 cm.

## 3 Conclusion

In this paper we studied the application of a magnetic levitation system to the seismic suspensions of long baseline GW antennas. From the results of the simulations shown in Sec. 2.2, we can summarize the following conclusions.

First, it appears that it is possible to levitate a mass comparable to the one of a Virgo-like payload, choosing reasonable dimensions for both the fixed and levitated magnet. This can also be done maintaining a free-gap of several centimeters, which would allow the insertion of a dielectric roof to separate the upper and the lower vacuum chambers.

A second interesting point is that, with a suitable magnet arrangement, the levitated system behaves also as a low frequency seismic attenuator, which improves thus the overall attenuation performance of the suspension. To improve the stability of the device, a pair of magnets aligned with opposing polarity can be used instead of a single one, as first suggested in Ref. . This configuration was validated experimentally in the Naples Virgo laboratory .

To give an example, we can consider as convenient configuration a fixed magnet of $`(40\times 20\times 10)\mathrm{cm}^3`$, a levitated one of $`(15\times 15\times 2)\mathrm{cm}^3`$ and a gap $`d=5`$ cm. In this configuration, the weights of the two magnets are 66.4 kg and 3.3 kg respectively, while the levitated payload $`M_p`$ is $`95`$ kg. As explained before, the horizontal oscillation frequency is, in this case, 0.23 Hz, giving an extra attenuation of about 74 db at 10 Hz.

The validity of the numerical model adopted for computing forces exerted between extended magnets has been experimentally tested with a small prototype suspension. The experimental results are in good agreement with the model. Of course, the study presented here only shows the feasibility of the principle; one of the main limitation is that we did not consider the coupling with the angular degrees of freedom of the suspended payload. A detailed study would be necessary to design a real suspension taking into account the need of controlling the other degrees of freedom of the mirror and other technical aspects, like the longitudinal control and the automatic alignment of the interferometer. In addition, a suspension designed for GW interferometers must not inject too many noise in the detector. Therefore, noise contributions originating from the magnetic configuration itself – like eddy and Johnson current effects – have to be taken into account on a full-scale prototype to see whether or not they limit the suspension performances. This is beyond the goal of this paper.