Import standard modules:
In [3]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
Out[3]:
Import section specific modules:
In [ ]:
HTML('../style/code_toggle.html')
So far, we've treated each individual measured visibility as a single measurement. Assuming a slowly and continuously rotating sky above the interferometer, we can turn the set of measured visibilities from different baselines into a set of samples of an underlying continuous visibility function $\mathcal{V}$. Under certain assumptions, we will see how this visibility function is related to the Fourier transform of the sky brightness distribution.
Based on what we have defined in $\S$ 4.3 ➞, 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$$
TLG:AC: Add image here which may help.
In the previous section, we have seen 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 (which are a function of the projected baseline) and modulated by the antenna response (assumed to be identical for all receivers so far):
$$\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, so that $\mathcal{V}_{\mathbf{b}}$ has dimensions of flux density:
Equation 4.4.1
In $\S$ 4.2 ➞, 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 ($u$,$v$,$w$)-coordinates, $\mathbf{s_0}$ defines the $w$ direction.
TLG:AC: Make a reference to chapter 3 (the direction cosines and to 4.2 prec where the $uvw$ coordinates are defined.
Therefore:
\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. 4.4.1 ⤵ can be expressed as a function of ($u$,$v$,$w$) and ($l$,$m$,$n$) as follows:
We will now express $d\Omega$ from Eq. 4.4.1 ⤵ in term of ($l$,$m$,$n$) coordinates. $d\Omega$ 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 using the Jacobian determinant, we end up with:
$$d \Omega = \frac{dl dm}{n} = \frac{dl dm}{\sqrt{1 - l^2 - m^2}}$$Equation 4.4.2
TLG:AC: Point to appendix for the proof.
If the following conditions are met:
i.e. $l,m << 1$ then:
$w(\sqrt{1-l^2-m^2}-1) \sim -\frac{1}{2}(l^2+m^2)w$
TLG:AC: This is too vague. The proof of the above is important. How the above conditions lead to the above eq + how does it help us to get the result below?
Eq. 4.4.2 ⤵ now 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$$
Equation 4.4.3
The Eq. 4.4.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 as:
$$ \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 in $\S$ 4.6 ➞
The Fourier relationship existing between $I_\nu$ and $\mathcal{V}$ is key, as it is the latter that is sampled by the interferometer. We remember that the interferometer is not sensitive to the sky but to the Fourier Transform 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 treat it as a continuous complex function.
$\mathcal{V}$ is a 2D function in 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 sample of the complex 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 Fourier Transform properties linking the visibility function to the brightness distribution, we will simply compute this integral numerically on a simple case. The Fourier Transform 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 Fourier basis 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.4.1 ⤵), through which the content of the sky is seen.
Figure 4.4.1: a single spatial frequency located at ($u$,$v$) and its associated fringe pattern on the sky.
For a 2-element interferometer, observing a sky which can be approximated as 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 with one another (Fig. 4.4.2 ⤵).
Figure 4.4.2: Relationship between the projected baseline, the ($u$,$v$,$w$) space and the ($l$,$m$,$n$) space.
From Eq. 4.4.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 from the observed field of view ($dldm$) coming from the direction $\boldsymbol{\sigma}(l,m)$. If the sky is filled with sources with complex structures, the visibility function being a function of ($u$,$v$,$w$), can be hard for a human to interpret.
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 an East-West baseline (see $\S$ 4.5.1 ➞).
We assume that the sky is only composed of a single extended source represented by a disk of unit brightness. For simplicity, this interferometer will observe this source at the phase center when the source is at transit.
We assume that this is the only source the interferometer sees in the sky, that the effect of the antenna pattern is negligible and that $w=0$.
The integral of Eq. 4.4.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 [8]:
# 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 [9]:
plotfringe(u=2,rad=0.2)
In Fig. 4.4.3 ⤵, the source is represented by the red disk over which we superimposed the fringe pattern. We see that an uneven fraction of the bright fringes cross the source. In this case, the resulting integral is positive. Let's try to increase the projected baseline size by increasing the value of $u$.
In [10]:
plotfringe(u=5,rad=0.2)
Figure 4.4.4: Same as in Fig. 4.4.3 ⤵ seen through a different fringe pattern with a longer projected baseline.
In Fig. 4.4.4 ⤵, the absolute value of the integral is $\sim$10 times lower than before, suggesting a more balanced contribution of dark and bright fringes to the integral (which is still 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 [11]:
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 [12]:
integral=plotintegral(UMAX=25.)
In Fig. 4.4.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:
$$ \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, is highly correlated with the geometry of the source we are observing.
Let's determine an approximate value for $u_\text{min}^{(1)}$.
In [13]:
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 [14]:
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}Indeed, in this example, the true object radius in unit of $l$ was 0.2. Our answer is thus a pretty good estimate!
In this simple example, we were able to measure the radius of an object. 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 that 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 resolution instruments.
See the following references for more details: Kraus, $\S$6-27 ⤴ and Taylor, $\S$ 16, p.338 ⤴.
Let's now use the visibility to determine the location of a source.
In order to reduce the bandwidth's impact on the fringe (recall that it causes the fringe pattern to be modulated by a sinc), the interferometer compensates for the delay $\tau$ by inserting a supplementary delay $\tau_c$ resulting in the shifting of the fringe pattern's phase center so that it coincides with the phase center. What happens if we remove do not implement any kind of fringe tracking and point the phase center towards the zenith? Well, we get a meridian transit interferometer.
For this example, we assume a single point source of unknown RA/DEC coordinates. We want to use the raw 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. This fringe will be modulated by the beam as the source crosses the array's field of view.
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.4.6 ⤵ in green) as the source passes through bright and dark fringes.
Figure 4.4.6: Derivation of the fringe rate and the transit time.
The Right Ascension can be determined if the transit time can be measured with high accuracy. The transit time can be derived from the maximum of the fringe enveloppe, which corresponds to the maximum elevation of the source multiplied by to the maximum antenna response of antennas pointing at the local meridian. TLG:RC: Please rewrite last sentence.EB:RC: I think that's what he meant? Not actually sure.. 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 rotation of the Earth.
TLG:GM: Check if the italic words are in the glossary.
In Fig. 4.4.6 ⤵, we estimate the fringe spacing to be $\sim$12 min ($\sim$ 5 periods are crossed in an 1h observation).
For a source located on the celestial equator ($\delta =0^\circ$), it will travel 15$^\circ$ per hour. This will correspond to a certain number of fringe spacings - counting them will thus give us information on the source's coordinates.
Since we know 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 estimate the declination of the source, we need to compare the fringe rate at the equator to the fringe rate of the fringe pattern in Fig. 4.4.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 known as fringe rate mapping (see Kogan, 1996 ⤴ and Walker, 1981 ⤴) exploits the information of the fringe rate offset compared to that at the phase center. TLG:GM: Check if the italic words are in the glossary. In a two-dimensional reference frame, the fringe rate offset can be defined 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.
By setting $m=f(l)$, Eq. 4.4.4 ⤵ becomes a linear function: $y=ax+b$
Each baseline produces a set of fringes on the sky with a different fringe spacing and orientation. This set of fringes rotates on the sky (due to 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. More distant sources will have a faster fringe rate contribution to the visibilities than 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 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 $m=f(l)$, giving all possible loci in the sky where the source can produce the same observed fringe rate.
If we have a second baseline, and thus a different measure of variation of the visibility as a function of time, we can derive a different ($\frac{du}{dt}$,$\frac{dv}{dt}$) parameter. This will give us another straight line.
The intersection of these lines will give the location of the source responsible for the variation of the two different baselines' measured amplitudes.
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 we have detected peaks, per baseline.
Ultimately, the collection of all straight line equations will produce intersections in the $(l,m)$ plane. These correspond to dominant sources in the field - without having to perform any imaging!
Having described the properties of the visibility function, we will address the sampling of the visibility function. To this end, we will focus on the variation of our measured ($u$,$v$,$w$) frequencies as a function of instrumental and observational parameters.
We have shown how a baseline can be expressed in a 3-D reference frame, and will next discuss how that baseline can be projected down to 2-D. Imagine being on the celestial sphere, looking down on to the baseline distribution. From this frame, the baselines look projected to a 2-D plane. Later in the chapter we will discuss the $\S$ 4.6 ➞ which allows us to construct an image of the sky from an approximate 2-D Fourier transform between the baseline visibility sampling and the sky. There are many caveats to this theorem, which will be discussed in good time. For now, the important thing to understand is: 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 the practice of radio interferometric observation, and will come up again when we discuss $\S$ 4.5.1 ➞.
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.