VO Science: The Relationship Between Clusters and their Galaxies (C. Miller 08/10/05)

In this lesson, we discuss how astronomers can integrate VO services into their day-to-day research activities.

A picture is worth a 1000 words:


Optical/X-ray image of the Coma Cluster

Coma Cluster in the Optical (Galaxies)

Coma Cluster in the X-ray (Gas)
www.damtp.cam.ac.uk/user/gr/public/gal_lss.html
 Credit: NASA/SAO/A.Vikhlinin et al.




 A movie is worth ?
Click here to see a movie of cluster formation (White et al.)


A High-z Cluster
Pieter van Dokkum (University of Groningen), ESA and NASA

What happens to the galaxies in clusters?

A nearby spiral galaxy

An old elliptical galaxy

Galaxy "collision"

http://quixote.uwaterloo.ca/~mbalogh/research/clusters.html
    
Click here to see a movie of a galaxy collision (Hibbard and Barnes)
Click here to see a movie of Ram Pressure Stripping (Vollner et al., Strasbourg)
“Strangulation”



Goto et al. 2003PASJ...55..757G

How do we quantify what we see? Add up the light (Broadly)





Bright, hot, young and blue stars

r-band (enhanced) older stars

Hα imaging (recent star-forming)







JHK Band Composite (old stars)

100μ Far Infrared (dust)

Radio (Syncotron from SN)

Bianchi et al., Condon et al. among others

How do we quantify what we see? Add up the light (Narrowly)


www.sdss.org

Line-Classification DiagramLine Classification Diagram

Miller et al. 2004

How do we quantify what we see? Measure the shapes of objects





Mostly Bulge (note the dusty disk)

Disk with a Bulge present

All Disk



What Figures do we publish?





Density-Morphology Relation

(Dressler et al. 1984)

Spectra-type vs. Density

(Miller et al. 2004)

Star-formation vs. Cluster Distance

Gomez et al. (2003)


How do we interpret what we've quantified?


Where do we get cluster data?  The Registry
Query on "galaxy clusters". Query on "clusters". Many Tabularskyservices (100s). 1 SIAP (CNOC), Few cone searches.
The Abell Cluster Catalog is arguably the most used cluster dataset.
The SDSS-C4 Galaxy Cluster Catalog is a recently published SDSS catalog using the spectroscopic sample.


Where do we get our galaxy data?
Cone Search Services (54 services with "galaxy" in the description)
Open SkyQuery (ADQL calls) (most nodes are galaxies)
Tabularskyservices

What can we do with the VO?

Putting it All Together: A specific Science Example

Below is an IDL script which utilizes most of the above VO calls. Additionally, it applies many IDL-specific astronomical libraries to the data found and returned by the VO services

Title:The Local Environment of Bright Cluster Galaxies (BCGs) and its Effect on their Properties
The local environment of a BCG provides insight into the (still not understood) formation process of the BCG itself. In terms of local environment, we will use the density of galaxies with in a small projected radius surrounding the BCG. We will also use the presence of any X-ray emission, both as a property of the local environment (if extended) or as a property of the BCG (if point-like). Infrared colors will be used as a mass proxy. Radio luminosity will be used as a BCG star-formation/AGN measure. The local environment will also include the cluster itself, and even the nearest cluster (as a function of cluster mass/luminosity). We'll want images (X-ray and optical) as well.

Plan: Using the C4 cluster Brightest Cluster Galaxy ID, find all galaxies to m_r = 23 and within a 150kpc radius.  Compare the shape, k-corrected luminosity, and k-corrected color of the BCG to the density of galaxies around the BCG. Look for any incidence of activity (star-formation and/or AGN) in the X-ray, and compare to the shape of the BCG and to the density of galaxies around the BCG. Measure the shape of the internal 100kpc of the clusters and compare to the 1-2Mpc shape of the clusters. Compare to the distance and direction of the nearest cluster. Use any radio data to classify the BCG as either star-forming or AGN. Use any available infrared colors as a proxy for the mass of the BCG (compared to models).

An inventory for the project

Functions/Tools:
          Luminosity Distances
          Absolute magnitudes
          K-corrections
          Angular size (on physical aperture)
          Statistical methods (means, KS-significance tests, etc)
          Data plotting routines
          Histograms
          Image display

Services:

SDSS DR2 SkyNode (or perhaps Cone Service)
ROSAT BSC and XSC source catalog
ROSAT Pointed catalog (WGACAT)
XMM SSC
2MASS XSC as either a SkyNode or Cone Search Service
FIRST radio matches (from SDSS or SkyNode)
SDSS SIAP
ROSAT SIAP
XMM SIAP

Data:

Cluster centers, masses/luminosties, shapes, with a BCG ID
SDSS 5 band galaxy magnitudes
Reddening
SDSS Shape measurements
X-ray fluxes, extents (ROSAT BSC and FSC)
X-ray fluxes (ROSAT PSPC)
X-ray fluxes (XMM)
2MASS J,H,K magnitudes
Radio (FIRST) flux          
Optical images (SDSS)
ROSAT All Sky Survey (RASS) Images
ROSAT PSPC Pointed images
XMM images


Let's do all of the above in 10 Steps.:


Step 1:  Declare variables, read in tables, etc.

Step 2: Get the SDSS data (catalogs and images).

Step 3: Get the SDSS-2MASS Cross-matches

Step 4: Do the k-corrections and compute the absolute magnitudes:

Step 5: Get the ROSAT All-Sky Survey data (catalogs and images)

Step 6: Get the ROSAT PSPC (hi-resolution) data (catalogs and images)

Step 7: Get the XMM data (catalogs and images)

Step 8: Save the BCG VO catalog data

Step 9: Do the science!

  1. Plot 

    plot, -(bcg_isophi[wkeep]-90),c4data[wkeep].ang1000, psym=4,xtitle="BCG Position Ange (degrees)", $
          ytitle="C4 Cluster Position Angle (within 1Mpc) (degrees)"

  2. Histogram

    h = histogram(abs(-(bcg_isophi[w]-90)-c4data[w].ang1000),omin=omin, binsize = 5)
    plot, 5*findgen(n_elements(h)) + omin, h, psym=10

  3. KS-test

    kstwo, delta_phi, zran,d,prob

  4. Spearman Correlation test

    result = r_correlate(-(bcg_isophi[w]-90), c4data[w].ang1000)

  5. Image viewers

    imdisp


The Code:
NOTE: You will need astro_lib (GSFC), kcorrect(Blanton), idlutils(SDSS) and separ.pro (see IDL Tutorial page) to make this work.
NOTE: You will need to copy and paste this file into a file called demo1.pro (demo1.pro).
NOTE: You will need to make a directory to store images called images/

PRO demo1, nosiap=nosiap
;If you don't want to download the images, type IDL> demo2,/nosiap 
  3. c4data = mrdfits('sdss_c4_dr2_published.fit',1,hdr)        ;Reads in cluster data
  4. radius_fixed = 1 50.0                                       ;Sets the aperture to 100kpc
  5. density = dblarr(n_elements(c4data.ra_mean))               ;Declare the density variable
  6. bcg_ur = dblarr(n_elements(c4data.ra_mean))                ;Declare the color variable
  7. bcg_absR = dblarr(n_elements(c4data.ra_mean))                ;Declare the absolute mag variable
  8. bcg_AB = dblarr(n_elements(c4data.ra_mean))                ;Declare the ellipticity variable
  9. bcg_isophi = dblarr(n_elements(c4data.ra_mean))            ;Declare the position angle variable
 10. bcg_absJ = dblarr(n_elements(c4data.ra_mean))                 ;Declare the J magnitude variable
 11. bcg_CPS = dblarr(n_elements(c4data.ra_mean))               ;Declare the RASS X-ray counts variable
 12. bcg_HR1 = dblarr(n_elements(c4data.ra_mean))               ;Declare the RASS X-ray hardness ratio
 13. bcg_PSPC_CTR = dblarr(n_elements(c4data.ra_mean))          ;Declare the PSPC X-ray counts variable
 14. bcg_EP = dblarr(n_elements(c4data.ra_mean))                ;Declare the XMM X-ray counts variable
;Loop over the cluster list
 16. FOR I = 0,200 DO BEGIN
 17.  radius =  zang(radius_fixed,c4data[I].z)/60.0   ;Convert 100kpc to angular size at redshift of the cluster
; Make the SIAP query to download the SDSS jpg images (url from idl> call_registry, 'SDSS', /SIAP, str=str)
 19.  IF not (keyword_set(nosiap)) THEN siapcall,c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, radius/60.0, $
 20.            url="http://casjobs.sdss.org/vo/DR4SIAP/SIAP.asmx/getSiapInfo?&FORMAT=image/jpeg&BANDPASS=*&", $
 21.            root="images/sdss_c4_"+strtrim(string(c4data[I].cluster_id),2)
; Specify the SkyPortal Client query 
 23.  qry = " SELECT o.ra,o.dec,o.expAB_r, o.isoPhi_r " + $
 24.        " FROM SDSSDR2:PhotoPrimary o " + $
 25.        " WHERE o.type = 3 AND o.petroMag_r < 23.0 " + $
 26.        " AND Region('Circle J2000 " + strtrim(string(c4data[I].ra_bcgphot,format='(f10.3)'),2) + $
 27.        " " + strtrim(string(c4data[I].dec_bcgphot, format='(f10.3)'),2) + $
 28.        " " +  strtrim(string(radius, format='(f4.2)'),2) + "') "
 29.  skyclient, qry=qry,str=sdss_gals                          ;Make the SkyPortal call
 30.  density[I] = double(n_elements(sdss_gals))/(!PI*(radius_fixed^2.0))   ;Calculate the surface density
;Find the BCG (as defined in the cluster data). The galaxy with the smallest separation.
 32.  separ, separ,  c4data[i].ra_bcgphot, sdss_gals.sdssdr2_ra,  c4data[i].dec_bcgphot, sdss_gals.sdssdr2_dec
 33.  min = min(separ, minit)
;Assign data to variables
 35.  bcg_AB[I] = sdss_gals[minit].sdssdr2_expAB_r;  Ellipticity
 36.  bcg_isoPhi[I] = sdss_gals[minit].sdssdr2_isophi_r; Position angle
 37.  chisq = 2.0
;Cross-match to 2MASS foor Infrared Magnitudes 
 39.  qry = " SELECT o.ra,o.dec, o.modelMag_u, o.modelMagErr_u, o.modelMag_g, o.modelMagErr_g, o.modelMag_r, " + $
 40.        " o.modelMagErr_r, o.modelMag_i, o.modelMagErr_i, o.modelMag_z,  o.modelMagErr_z, " + $
 41.        " o.extinction_u, o.extinction_g, o.extinction_r, o.extinction_i, o.extinction_z, " + $
 42.        " t.j_m, t.k_m, t.h_m, t.j_msigcom, t.k_msigcom, t.h_msigcom " + $
 43.        " FROM SDSSDR2:PhotoPrimary o, TWOMASS:PhotoPrimary t " + $
 44.        " WHERE XMATCH(o,t)<" + strtrim(string(chisq),2) + " " + " AND o.type = 3 " + $
 45.        " AND Region('Circle J2000 " + strtrim(string(c4data[I].ra_bcgphot,format='(f10.3)'),2) + $
 46.        " " + strtrim(string(c4data[I].dec_bcgphot, format='(f10.3)'),2) + $
 47.        " " +  strtrim(string(radius, format='(f4.2)'),2) + "') "
 48.  skyclient,qry=qry,str=sdss_2mass_gals
 49.  separ, separ,  c4data[i].ra_bcgphot, sdss_2mass_gals.sdssdr2_ra, $
 50.                 c4data[i].dec_bcgphot, sdss_2mass_gals.sdssdr2_dec
 51.  min = min(separ, minit)
 52.  bcg_absR[I] = sdss_2mass_gals[minit].sdssdr2_modelMag_r - 5*alog10(lumdist(c4data[I].z,/silent)*1e6) + 5
 53.  bcg_absJ[I] = sdss_2mass_gals[minit].twomass_j_m - 5*alog10(lumdist(c4data[I].z,/silent)*1e6) + 5
;If you have kcorrect installed, kcorrect the magnitudes:
 55.  mags = dblarr(8,n_elements(sdss_2mass_gals))
 56.  magerrs = dblarr(8,n_elements(sdss_2mass_gals))
 57  kcorrct = dblarr(8, n_elements(sdss_2mass_gals))
 58.  FOR K = 0, n_elements(sdss_2mass_gals) -1 DO BEGIN
 59.   mags[*,K] = [sdss_2mass_gals[K].sdssdr2_modelMag_u - sdss_2mass_gals[K].sdssdr2_extinction_u, $
 60.                sdss_2mass_gals[K].sdssdr2_modelMag_g - sdss_2mass_gals[K].sdssdr2_extinction_g, $
 61.                sdss_2mass_gals[K].sdssdr2_modelMag_r - sdss_2mass_gals[K].sdssdr2_extinction_r, $
 62.                sdss_2mass_gals[K].sdssdr2_modelMag_i - sdss_2mass_gals[K].sdssdr2_extinction_i, $
 63.                sdss_2mass_gals[K].sdssdr2_modelMag_z - sdss_2mass_gals[K].sdssdr2_extinction_z, $
 64.                sdss_2mass_gals[K].twomass_j_m, sdss_2mass_gals[K].twomass_h_m, sdss_2mass_gals[K].twomass_k_m]
 65.   magerrs[*,K] = [sdss_2mass_gals[K].sdssdr2_modelMagErr_u, sdss_2mass_gals[K].sdssdr2_modelMagErr_g, $
 66.                  sdss_2mass_gals[K].sdssdr2_modelMagErr_r, sdss_2mass_gals[K].sdssdr2_modelMagErr_i, $
 67.                  sdss_2mass_gals[K].sdssdr2_modelMagErr_z, sdss_2mass_gals[K].twomass_j_msigcom, $
 68.                 sdss_2mass_gals[K].twomass_h_msigcom, sdss_2mass_gals[K].twomass_k_msigcom]
 69.  ENDFOR
 70.  filterlist = ['sdss_u0.par', 'sdss_g0.par', 'sdss_r0.par', 'sdss_i0.par', 'sdss_z0.par', $
 71.                'twomass_J.par', 'twomass_H.par', 'twomass_Ks.par']
 72.  kcorrect,mags, magerrs, dblarr(n_elements(sdss_2mass_gals)) + c4data[I].z , kcorr, filterlist=filterlist, band_shift=0.0
 73.  bcg_absR[I] = bcg_absR[I] - kcorr[2,minit]
 74.  bcg_absJ[I] = bcg_absJ[I] - kcorr[5,minit]
 75.  bcg_ur[I] = (mags[0,minit] - kcorr[0,minit]) - (mags[2,minit] - kcorr[2,minit])
;Cross-match to ROSAT BSC/FSC 
 77.  qry = " SELECT o.ra, o.dec, o.srcCps, o.srcCpsErr, o.ext, o.hr1, o.hr2 " + $
 78.        " FROM Rosat:photoprimary o " + $
 79.        " WHERE Region('Circle J2000 " + strtrim(string(c4data[I].ra_bcgphot,format='(f10.3)'),2) + $
 80.        " " + strtrim(string(c4data[I].dec_bcgphot, format='(f10.3)'),2) + $
 81.        " " +  strtrim(string(radius, format='(f4.2)'),2) + "') "
 82.  skyclient,qry=qry,str=rosat_gals
 83.  separ, separ,  c4data[i].ra_bcgphot, rosat_gals.rosat_ra,  c4data[i].dec_bcgphot, rosat_gals.rosat_dec
 84.  min = min(separ, minit)
 85.  bcg_CPS[I] = rosat_gals[minit].rosat_srcCps
 86.  bcg_HR1[I] = rosat_gals[minit].rosat_hr1
;If there is a source, get the image
 88.  IF (bcg_CPS[I] ne 0) THEN BEGIN
 89.   IF not (keyword_set(nosiap)) THEN siapcall,c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, 0.1, $
 90.         url="http://skyview.gsfc.nasa.gov/cgi-bin/vo/sia.pl?survey=rass-cnt&",$
 91.         root="images/rosat_bcg_rass_c4_"+strtrim(string(c4data[I].cluster_id),2);,/metadata,str=str
 92.  ENDIF
;Call the WGACAT Cone Search (discovered via idl> call_registry, 'wgacat', str=str)
 94.  conecall, c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, radius/60.0, str=str, $
 95.            url = "http://heasarc.gsfc.nasa.gov/cgi-bin/vo/cone/coneGet.pl?table=wgacat&r"
 96.  sizeit = size(str)
;If there is a WGACAT pointed observation, get the image
 98.  IF (sizeit[1]  gt 1) THEN BEGIN
 99.   separ, separ,  c4data[i].ra_bcgphot, str.ra,  c4data[i].dec_bcgphot, str.dec
100.   min = min(separ, minit)
101.   bcg_PSPC_CTR[I] = str[minit].count_rate
102.   IF not (keyword_set(nosiap)) THEN siapcall,c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, 0.1, $
103.           url="http://skyview.gsfc.nasa.gov/cgi-bin/vo/sia.pl?survey=%5epspc&", $
104.           root="images/rosat_pspc_c4_"+strtrim(string(c4data[I].cluster_id),2);,/metadata, str=str
105.   IF not (keyword_set(nosiap)) THEN siapcall,c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, 0.1, $
106.           url="http://skyview.gsfc.nasa.gov/cgi-bin/vo/sia.pl?survey=hri&", $
107.           root="images/rosat_bcg_hri_c4_"+strtrim(string(c4data[I].cluster_id),2);,/metadata,str=str
108   ENDIF
;Call the XMM Serendipitous Source Catalog Cone Search
110.  conecall, c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, radius/60.0, str=str, $
111.            url = "http://heasarc.gsfc.nasa.gov/cgi-bin/vo/cone/coneGet.pl?table=xmmssc&"
112.  sizeit = size(str)
;If there is XMM-SSC data, get the image
114.  IF (sizeit[1] gt 1) THEN BEGIN
115.   separ, separ,  c4data[i].ra_bcgphot, str.ra,  c4data[i].dec_bcgphot, str.dec
116.   min = min(separ, minit)
117.   bcg_EP[I] = str[minit].ep_flux
118.   IF not (keyword_set(nosiap)) THEN  siapcall,c4data[I].ra_bcgphot, c4data[I].dec_bcgphot, 0.1, $
119.           url="http://xsa.vilspa.esa.es:8080/aio/jsp/siap.jsp", $
120.           root="images/xmm_bcg_c4_"+strtrim(string(c4data[I].cluster_id),2);, /metadata,str=str
121.   ENDIF
122. ENDFOR
;Write everything out to a file 
124. openw,1,'bcg_data.dat'
;    FOR I = 0, n_elements(density)- 1 DO BEGIN
126  FOR I = 0, 200 DO BEGIN 
127.    printf,1,density[I], bcg_ur[I], bcg_AB[I], bcg_isophi[I],bcg_absR[I], $
128.             bcg_absJ[I], bcg_CPS[I], bcg_EP[I], bcg_PSPC_CTR[I], $
129.             format = '(9(1x, f10.5))'
130. ENDFOR 
131. close, 1
132. stop
133. END


To Do the Science we must make our figures and run our statistics:

PRO demo1_science
rdfloat, 'bcg_data.dat', density, bcg_ur, bcg_AB, bcg_isophi, bcg_absR, $
         bcg_absJ, bcg_CPS, bcg_EP, bcg_PSPC_CTR,/double
c4data = mrdfits('sdss_c4_dr2_published.fit',1,hdr)        ;Reads in cluster data
;First, "fix" the C4 Cluster position angles:
w = where(c4data.ang1000 gt 90)
c4data[w].ang1000 = c4data[w].ang1000 -180
;Next, select only those clusters and BCGs that are in fact elliptical and bright
;Choose only clusters above z = 0.05 to avoid elliptical selection effects.
wkeep = where(bcg_AB lt 0.65 and c4data.semiminor1000/c4data.semimajor1000 lt 0.65 and bcg_AB gt 0.0 and $
        bcg_absR lt -22.0 and c4data.z ge 0.05)                             
sorted = sort(density[0:200])
!P.MULTI=[0,2,2]
set_plot, 'ps'
device,/inches,ysize=7.0,scale_factor=1.0
device,/inches,xsize=8.0,scale_factor=1.0
device, filename="bcg_v_cluster_posang.ps"
plot, -(bcg_isophi[0:200]-90),c4data[0:200].ang1000, psym=4,xtitle="BCG P.A. (degrees)", $
      ytitle="C4 Cluster P.A. (within 1Mpc) (degrees)", xr=[-90,90], yr = [-90,90]
oplot, -(bcg_isophi[sorted[190:200]]-90),c4data[sorted[190:200]].ang1000, psym=2  
;correlate
delta_phi =  abs(-(bcg_isophi[wkeep]-90)-c4data[wkeep].ang1000)
wtmp = where(delta_phi gt 90)
delta_phi[wtmp] = 180 - delta_phi[wtmp]
zran = dblarr(100, n_elements(wkeep))
probd = dblarr(100)
FOR I = 0, 99 DO BEGIN zran[I,*] = randomu(iseed, n_elements(wkeep))*90.0
FOR I = 0, 99 DO BEGIN kstwo, delta_phi, zran[I,*],d,prob & probd[I] = prob
print, "The probability that the difference between the BCG and the Clustet PAs is random: ", median(probd)
plot, density[0:200], bcg_AB[0:200], psym=4,$
      xtitle="Local BCG density", ytitle="BCG Ellipticity"
oplot,density[sorted[190:200]], bcg_AB[sorted[190:200]], psym=2
plot, density[0:200], bcg_ur[0:200], psym=4,yr=[2,4], $
      xtitle='Local BCG density', ytitle="BCG color"
oplot,density[sorted[190:200]], bcg_ur[sorted[190:200]], psym=2  
plot, density[0:200], bcg_absJ[0:200], psym=4,yr=[-24,-21], $
      xtitle="Local BCG Density", ytitle="BCG J-band Luminosity"
oplot,density[sorted[190:200]], bcg_absJ[sorted[190:200]], psym=2 
device,/close
!P.MULTI=[0,1,1]
set_plot, 'ps'
device,/inches,ysize=7.0,scale_factor=1.0
device,/inches,xsize=8.0,scale_factor=1.0
device, filename="bcg_xray_histo.ps"
w = where(bcg_CTR gt 0)
plothist, density, bin=0.0001                   
plothist, density[w], bin=0.0001,/overplot,/fill
device,/close
;You'll need the imdisp routine for this to work:
!P.MULTI=[0,3,4]
w = where(bcg_PSPC_CTR gt 0)
set_plot, 'ps'
device,/inches,ysize=7.0,scale_factor=1.0
device,/inches,xsize=8.0,scale_factor=1.0
device, filename="bcg_images.ps"
FOR I = 0, 3 DO BEGIN
  name1 = 'images/rosat_pspc_c4_' + strtrim(string(c4data[w[i]].cluster_id),2) + '_image1.gif'
  name2 = 'images/sdss_c4_' + strtrim(string(c4data[w[i]].cluster_id),2) + '_image0.jpg' 
  name3 = 'images/rosat_bcg_rass_c4_' + strtrim(string(c4data[w[i]].cluster_id),2) + '_image1.gif' 
  read_gif,name1,image1
  read_jpeg,name2,image2
  read_gif,name3,image3
  imdisp, image1
  imdisp, image2
  imdisp, image3
ENDFOR
stop
END