In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
In [ ]:
HTML('../style/code_toggle.html')
Assuming a slowly and continuously rotating sky above the interferometer, the set of measured visibilities from different baselines is a set of samples of the continuous visibility function $\mathcal{V}$. Under certain assumptions, we will see how this visibility function is related to the Fourier transform (FT) of the sky brightness distribution.
Based on what we have defined in $\S$ 4.2 ➞, we will focus on the visibility function which is a caracteristic of the sky, as seen through an interferometer. It is related to the contrast in the fringe pattern created by the combination of signals received by the antennas on each baseline.
As seen in $\S$ 1.2 ➞, the specific intensity $I_\nu$ (or the surface brightness) is defined as the received power per unit solid angle, per unit frequency interval, per unit collecting area: $$ dP_{\nu} = I_{\nu} d\Omega d\nu d A_{\text{eff}} $$ where $dP_{\nu}$ is in units of W and $I_{\nu}$ is in units of W m$^{-2}$ sr$^{-1}$ Hz$^{-1}$.
In the integrated form, the received power $P_{\text{rec}}$ from a source with flux density $S$ over the bandwidth $\Delta \nu$, flowing through a collecting area $A$ is: $$P_{\text{rec}}=\frac{1}{2} A S \Delta \nu$$
At a location $\boldsymbol{\sigma}$ from the phase center $\mathbf{s_0}$, the region of the sky defined by $d\Omega$ will contribute to the received power over bandwidth $\Delta \nu$ and collecting surface $A_{\text{eff}}$ as: $$dP(\boldsymbol{\sigma})= \frac{1}{2} A_{\text{eff}}(\boldsymbol{\sigma})I_\nu(\boldsymbol{\sigma})\Delta\nu d\Omega$$
We have seen in the previous section that the quantity measured by an interferometer is the sum of the contributions from different parts of the sky within the solid angle $\Omega$. The complex visibility includes the measurements of both the odd part and the even part of the specific intensity $I_\nu$, as seen through spatial filters (depending on the projected baseline) and modulated by the antenna response (we assumed so far that these responses were identical for all receivers):
$$\boxed{\boxed{\mathcal{ V_{\mathbf{b}} } = \int_{\Omega} A_\text{eff}(\boldsymbol{\sigma})I_\nu(\boldsymbol{\sigma}) e^{-\imath 2\pi \frac{\textbf{b}\cdot \boldsymbol{\sigma}}{\lambda}} d\Omega}}$$We will use an expression of the visibility where the effective area has been normalized to give $\mathcal{V}_{\mathbf{b}}$ the dimension of a flux density:
In $\S$ 4.1 ➞, we defined various coordinates systems to represent the baseline in a sky-friendly reference frame. We will now use them to have an explicit expression of the visibility $\mathcal{V}_{\mathbf{b}}$.
$\boldsymbol{\sigma}$ is the direction difference vector defined as $\boldsymbol{\sigma}= \mathbf{s} - \mathbf{s_0}$.
In the ($u$,$v$,$w$) coordinate system, $\mathbf{s_0}$ define the $w$ direction. We therefore have :
\begin{eqnarray} \mathbf{s_0}&=& \begin{pmatrix} 0 \\ 0 \\ 1 \\ \end{pmatrix} \end{eqnarray}\begin{eqnarray} \mathbf{\boldsymbol{\sigma}} &=& \begin{pmatrix} l \\ m \\ n \\ \end{pmatrix} \end{eqnarray}\begin{eqnarray} \mathbf{b}_{\lambda} &=& \mathbf{ \frac{\mathbf{b}}{\lambda}} = \begin{pmatrix} u \\ v \\ w \\ \end{pmatrix} \end{eqnarray}The scalar product $\mathbf{b} \cdot \boldsymbol{\sigma}$ in Eq. 1 ⤵ can therefore be expressed as a function of ($u$,$v$,$w$) and ($l$,$m$,$n$) as following:
We will now express the quantity $d\Omega$ in Eq. 1 ⤵ in term of ($l$,$m$,$n$) coordinates. This quantity describes an element of solid angle on the celestial sphere. As $d \Omega = d\theta \sin \theta d\phi$, we may also express the surface element in terms of ($l$,$m$,$n$). By computing the Jacobian determinant, we end up with:
$$d \Omega = \frac{dl dm}{n} = \frac{dl dm}{\sqrt{1 - l^2 - m^2}}$$If the following conditions are met:
i.e. $l,m << 1$ then we can take simplify the expression of $\mathcal{V}$ by keeping the first order development of the argument of the exponential:
$w(\sqrt{1-l^2-m^2}-1) \sim -\frac{1}{2}(l^2+m^2) \sim \mathcal{O}(l,m) $
The integral becomes:
$$ \mathcal{V}(u,v,w \sim 0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{A(l,m) I_\nu (l,m)}{\sqrt{1 - l^2 - m^2}} e^{ -\imath 2\pi (ul+vm)}dl dm$$
The Eq. 3 ⤵ is no longer a function of $w$, and it now takes the shape of a 2-D Fourier transform (see $\S$ 2.4 ➞) with (u,v) being the Fourier pairs of (l,m).
The inverse transform can be written:
$$ \frac{A(l,m) I_\nu(l,m)}{\sqrt{1 - l^2 - m^2}} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathcal{V}(u,v) e^{ +\imath 2\pi (ul+vm)}du dv$$$$ \mathcal{V}(u,v) \sim \mathscr{F} \{I_\nu \}(u,v) $$We will come back to the consequences of this form later, in $\S$ 4.5 ➞
The Fourier relationship existing between $I_\nu$ and $\mathcal{V}$ if of primary interest as the latter is sampled by the interferometer. We remember that the interferometer is not sensitive to the sky but to the FT of the sky. To recover some knowledge of the sky through $\mathcal{V}$, one needs to understand the general shape of the visibility function. We will consider it as a continuous complex function.
$\mathcal{V}$ is a 2D function lying in the Fourier space. At a time $t$ and frequency $\nu$, one baseline and one direction $\mathbf{s_0}$, will provide a unique ($u$,$v$) point in the fourier space, which corresponds to one complex sample of the visibility function.
In the earlier days of interferometry, only few samples where available and the inspection of the samples in a 2D plane was unpractical. We usually inspect the amplitude of the visibility samples as a function of the uv distance $r_{uv}$ defined:
$$ r_{uv} = \sqrt{u^2+v^2}$$Instead of exploiting the FT properties linking the visibility function to the brightness distribution, we will simply compute this integral numerically on a simple case. The FT will be more extensively used in $\S$ 5.1 ➞ for imaging.
The visibility is expressed as the integral of the product of the intensity distribution with a complex exponential which is called the Fourier kernel.
Using a vector definition of the intensity distribution as a function of direction, this integral can be interpreted as the scalar product between the intensity distribution vector $\mathbf{I_\nu}$, and the Fourier basis kernel function $\mathbf{f}_{u,v}^{l,m}$:
$$ \mathcal{V}= \langle \mathbf{I_\nu} \cdot \mathbf{f}_{u,v}^{l,m}\rangle$$with $\mathbf{f}_{u,v}^{l,m}= e^{-2j\pi (ul+vm)}$
This can be seen as the projection of $\mathbf{I_\nu}$ on the basis vector $\mathbf{f}_{u,v}^{l,m}$. The complex visibility of the baseline with ($u$,$v$) coordinates, is therefore the coefficient of the intensity distribution projected on the basis Fourier vector of frequency ($u$,$v$). This operation therefore filters the magnitude of spatial frequency ($u$,$v$) contained in the intensity distribution.
The Fourier vector can be seen as a fringe pattern projected on the sky (Fig. 4.3.1 ⤵), through which the content of the sky is seen.
Figure 4.3.1: a single spatial frequency located at ($u$,$v$) and its associated fringe pattern on the sky.
Considering the case of a 2-element interferometer, observing a sky which can be locally be approximated by a plane, we can link the physical baseline $\mathbf{b}$, the projected baseline $\mathbf{b_\text{proj}}$ in the ($u$, $v$, $w$) frame, the phase center $\mathbf{s_0}$ and a direction of observation $\mathbf{\sigma}$ in the ($l$, $m$, $n$) frame (Fig. 4.3.2 ⤵).
Figure 4.3.2: Relationship between the projected baseline, the ($u$,$v$,$w$) space and the ($l$,$m$,$n$) space.
From Eq. 3 ⤵, for a single projected baseline with coordinates ($u$,$v$,$w$), the value of the visibility function $\mathcal{V}(u,v,w)$ is the sum of all contributions observed the field of view ($dldm$) coming from the direction $\boldsymbol{\sigma}(l,m)$. If the sky is filled with sources with complex structures, the shape of the visibility function as a function of ($u$,$v$,$w$) can be hard to interpret and some modelling is required.
Consider a 2-element interferometer projecting a baseline ($u$,$v=0$,$w=0$) associated with a set of fringes along the $m$ axis (for simplification). This is possible with a East-West baseline (see $\S$ 4.4 ➞).
We assume that the sky is only composed of a single extended source represented by a disk of unit brightness. For simplificity, this interferometer will observe this source at the phase center and during the transit.
We assume that the sky is only composed of this source, that the effect of the antenna pattern is negligible and that $w=0$.
The integral of Eq. 3 ⤵ therefore reduces to computing the integral of the fringe pattern over the source: $$ \mathcal{V} = \int_{\text{disk}} e^{-2j\pi ul}dl$$
In [ ]:
# 1 East-West baseline observing a disk at the phase center
from matplotlib.patches import Circle
def plotfringe(u=4,rad=0.2):
global radius
radius=rad
# preparing (l,m,n) space
Npointsl=1001
ll=np.linspace(-1.,1.,Npointsl)
l,m=np.meshgrid(ll,ll)
# Definition of the disk
#radius=0.1234 # angular radius of the object in l,m coordinates
# projected baseline length on the u axis
#u=4
# generate fringe pattern
tabcos=np.real(np.exp(-2j*np.pi*u*l))
# plotting the fringe pattern and the source
pxrad=radius*Npointsl/2
circle=Circle((500,500),pxrad,color='r',alpha=0.5,fill=True)
fig,ax =plt.subplots(figsize=(6,6))
im=plt.imshow(np.abs(tabcos),interpolation=None,cmap="winter")
ax.add_patch(circle)
#center=l[(Npoints-1)/2,(Npoints-1)/2]
# Compute the absolute value of the integral of the fringe over the source
w=np.where(np.sqrt(l**2+m**2) <= radius)
integral=np.sum(tabcos[w])
print "Integral="+str(integral)
In [ ]:
plotfringe(u=2,rad=0.2)
In Fig. 4.3.3 ⤵, the source is represented with the red disk on top of which the fringe pattern is super imposed. We see that an uneven fraction of the bright fringes is crossing the source. The resulting integral is positive in this case. Let's try to increase the projected baseline size by increasing the value of $u$.
In [ ]:
plotfringe(u=5,rad=0.2)
Figure 4.3.4: Same as inf Fig. 4.3.3 ⤵ seen through a different fringe pattern with a longer projected baseline.
In Fig. 4.3.4 ⤵, the absolute value of the integral is $\sim$10 times lower than before, suggesting a more balanced contribution of dark and white fringes to the integral (the latter being slightly dominated by the dark fringes).
We can understand that, as the width of the fringes decreases, the integral will ultimately statistically converge towards 0. Indeed the probability of evenly covering the source with the same fraction of bright and dark fringes increases.
Let's focus on the variation of the absolute value of the integral, as a function of increasing $u$.
In [ ]:
def plotintegral(UMAX=15):
%matplotlib inline
global u
from matplotlib.patches import Circle
#UMAX=5. # adjust it to larger values if no zeroes is encountered in next plot
Npointsl=1001
Npointsu=500
ll=np.linspace(-1.,1.,Npointsl)
l,m=np.meshgrid(ll,ll)
u=np.arange(Npointsu)*UMAX*1./Npointsu
w=np.where(np.sqrt(l**2+m**2) <= radius)
integral=np.array([])
for du in u:
tabcos=np.real(np.exp(2j*np.pi*du*l))
integral=np.append(integral,np.abs(np.sum(tabcos[w])))
normintegral=integral/np.max(integral)
plt.xlabel('Spatial frequency u')
plt.ylabel('Normalized integral over source')
plt.plot(u,normintegral,".-")
return normintegral
In [ ]:
integral=plotintegral(UMAX=25.)
In Fig. 4.3.5 ⤵, from $u=0$ to $u=25$, we can notice that the integral is close to zero at specific values of $u$ (denoted $u_\text{min}^{(n)}$). This corresponds to the fringe spacing where the integral of the source over the dark fringes is equal to the integral over the bright fringes.
Namely, when
$$ \int_{\text{bright }\cap\text{ source}} e^{-2j\pi ul}dl \approx \int_{\text{dark }\cap\text{ source}} e^{-2j\pi ul}dl$$A null in the integral (i.e. of the visibility function) corresponds to the case where the contrast of the fringe (over the source) is zero. In this particular case, we say that the source is resolved. The first value $u_\text{min}^{(1)}$, where the integral is minimum, has therefore a strong relationship with the geometry of the source we are observing.
Let's determine an approximate value $u_\text{min}^{(1)}$.
In [ ]:
def findumin(normintegral,ulim=5):############
# Adjust ulim to search for the first minimum of the integral
#
#ulim=5 # should be an value larger than the first minimum
wloc=np.where(u <= ulim)
############
locmin=np.where(normintegral[0:wloc[0][-1]] == np.min(normintegral[0:wloc[0][-1]]))
print "Index first minimum = "+str(locmin[0][0])
print "Normalized integral value at first min = "+str(normintegral[locmin][0])
print "Spatial frequency at first min = "+str(u[locmin][0])
umin=u[locmin][0]
deltal=1.22/(2*umin) # Bessel function : J1(3.83)=0=J1(2*pi*f*r)
print
print "Spatial scale at first min = "+str(deltal)
print "True object radius = "+str(radius)
In [ ]:
findumin(integral,ulim=5)
We identified the value of $u_\text{min}^{1}$ to be 3.05. As the object is a disk, its visibility function will be described by the first order Bessel function $\mathcal{J}_1$. We know that $\mathcal{J}_1(2\pi u_\text{min}^{1} l)=\mathcal{J}_1(3.8317) \approx 0$. We encounter a first null when:
\begin{eqnarray} 2 \pi u l &\approx& 3.8317 \\ l&\approx&\frac{1.22}{2 u_\text{min}^{1}} \\ l&\approx&0.2 \end{eqnarray}As we set in this example, the true object radius in unit of $l$ was 0.2.
With this small example, we performed the measurement of the radius of an object along the projected baseline. If you were observing a remote star with a sufficiently good system and a relatively long baseline at a short wavelength, it would be possible to measure the radius of this star, without any imaging processing.
This example illustrates how interferometrists were (and still are!) able to resolve far-away objects with high angular resolution by combining the signals from low-angular resolutions instruments.
Let's move on to use of the visibility to determine the location of the source.
For more study of the shape of the source from the visibility, please refer to Kraus, $\S$6-27 ⤴ and Taylor, $\S$ 16, p.338 ⤴.
In order to reduce the effect of the bandwidth on the fringe (i.e. causing the fringe pattern to be tampered by an sinc), the interferometer compensates for the delay $\tau$ by inserting a supplementary delay $\tau_c$ resulting in the shifting of center of the fringe pattern in the sky to follow the direction of the phase center. Let's remove this tracking for the moment and let's assume that the phase center is directed towards the zenith. We have a meridian transit interferometer.
For this example, we assume a single point source of unknown RA/DEC coordinates. We want to use the measurement of the interferometer to accurately determine its coordinates. We also assume that the field of view of the antenna elements is sufficiently large to vary smoothly accross the observation. We point the antennas toward the local meridian ($H=0^\text{h}$). The transit interferometer will project a fringe pattern on the sky which will be modulated by the beam pattern as the source cross the array factor and the element pattern.
A source on the celestial sphere will cross the projected fringe pattern and will create a variation in the interferometer response as a function of the LST (Fig. 4.3.6 ⤵ in green) as the source will pass in bright and dark fringes.
Figure 4.3.6: Derivation of the fringe rate and the transit time.
The Right Ascension can be determined if the transit time can be measured with a high precision. The transit time can be derived from the maximum of the fringe enveloppe, corresponding to the maximum elevation of the source combined to the maximum of the response of the antennas pointing at the local meridian. The transit time correspond to the moment when the LST is equal to the RA (according to $\S$ 3.2 ➞). In our example, the RA is $\sim$$13^\text{h}07^\text{m}$.
To estimate the declination, we need to measure the fringe spacing and the fringe rate of the source using the Earth Rotation.
On Fig. 4.3.6 ⤵, we estimate the fringe spacing to be $\sim$12 min ($\sim$ 5 periods are crossed in 1h of observation).
For a source located on the celestial equator ($\delta =0^\circ$), it will travel 15$^\circ$ per hour and a certain number of fringe spacings that we need to compute. Knowing the fringe spacing, we can derive the fringe rate at the equator with $\frac{d\phi}{d\theta}\Bigr|_\text{eq}=\frac{\Delta l_f}{15^\circ \text{per h}} \approx 0.77 h \approx 4^\text{m}33^\text{s}$.
Conversely, a source close to the NCP ($\delta \sim 90^\circ$) will have a very small fringe spacing and fringe rate.
To have an estimation of the source declination, we need to compare the fringe rate at the equator to that measured on the fringe pattern of Fig. 4.3.6 ⤵.
$\cos \delta= \frac{d\phi}{d\theta}\Bigr|_\text{eq} / \frac{d\phi}{d\theta}\Bigr|_\text{mes}=\frac{4^\text{m}33^\text{s}}{12^\text{m}}=0.3825 \leftrightarrow \delta \approx 67.7^\circ$.
A technique, called fringe rate mapping (see Kogan, 1996 ⤴ and Walker, 1981 ⤴) exploits the information of the fringe rate offset compared to that at the phase center. In a two-dimensional reference frame, the fringe rate offset can be written as:
$$ \omega_{\text{frm}}= 2\pi (\frac{du}{dt}l+\frac{dv}{dt}m)$$
where ($u$,$v$) are the coordinates of the projected baseline and ($l$,$m$) the direction cosine coordinates from the phase center.
Eq. 4 ⤵ is an equation of a straight line $y=ax+b$ once rearranged as $m=f(l)$.
Each baseline produces a set of fringes on the sky with different fringe spacing and orientation. This set of fringes rotates on the sky (as a result of the rotation of the baseline as seen from the source) and each source will modulate the visibility depending on its distance to the phase center. Far away sources will have a faster fringe rate contribution to the visibilities that sources close to the phase center.
If the sky is composed of only one source at an unknown position, the visibility amplitude will have a periodic behavior at a given the fringe rate (see previous example).
If we can measure the periodic behavior in the visibility plane as a function of time (i.e. measuring $\frac{du}{dt}$ and $\frac{dv}{dt}$), then we can draw the straight line of equation$m=f(l)$, giving all possible locii in the sky where the source can produce the same observed fringe rate.
If we have a second baseline, giving a different measure of variation of the visibility as a function of time, we will derive different ($\frac{du}{dt}$,$\frac{dv}{dt}$) parameter which will give another equation of straight line.
The intersection of the two lines (equivalent to solving the one equation with the two unknowns $l$,$m$) will give the locus of the source responsible for the variation on the two different baseline.
If we have multiple bright sources, one must study the spectrum of the variation of the visibility and search for the various distribution of peaks (each giving values of $\frac{du}{dt}$ and$\frac{dv}{dt}$). By taking the FFT of each visibility as a function of time we can derive as many lines as detected peak, per baseline.
Ultimately, the collection of all straight line equations will produce intersection in the l,m plane, corresponding to the dominant sources without having to image the source.
Now that the properties of the visibility function have been described, we will address the sampling of the visibility function by focusing on the variation of the measured ($u$,$v$,$w$) frequency, as a function of instrumental and observational parameters.
We have shown how a baseline can be transformed into a 3-D reference frame, we will introduce how that baseline is projected down to 2-D. We can imagine being at a point on the celestial sphere looking down on to the baseline distribution. From this point, the baselines look projected to a 2-D plane. Later in the chapter we will discuss the Van Cittert-Zernike Theorem ➞ which allows an image of the sky to be constructed from an approximate 2-D Fourier transform between the baseline visibility sampling and the sky. There are many caveats to this theorem, but for now we are getting ahead of ourselves. At the moment what is important to understand is that for a given point on the celestial sphere the 3-D baselines can be projected onto a 2-D plane, and as the Earth rotates (and the Celestial sphere stays fixed) this changes the projection of each baseline. This change in projection is key to how we observe in radio astronomy, and will come up again later when we discuss uv coverage ➞.
Important things to remember
• The visibility function reduces to the 2D Fourier transform of the intensity distribution under simplyfing assuptions.
• The visibility integral at ($u$,$v$) can be seen as an operator which returns the coefficient of the spatial frequency ($u$,$v$).
• By studying the properties of this integral, one can derive the angular size and position of a source without imaging.