When multiple screens are present, the number of paths from source to observer
become much larger. Below, we derive the paths taken for linear screens. The
result is implemented in the Screen1D class, and
this class is used in the tutorial on Simulating a multiple-screen system. To
help visualize how the resulting wavefield depends on the properties of the
screen, have a look at examples/two_screen_interaction.py.
Multiple screens
Consider screens with linear features between the telescope and the
pulsar. Consider a coordinate system in which z is along the line of
sight (with direction \(\hat{z}\)), and the lines on a given screen \(s\)
describing the features are given by,
(1)\[ d_{s}\hat{z} + p \hat{r} + s \hat{u}\]
where \(p \hat{r}\) is a cylindrical radius from the line of sight to the line
(ie., \(\hat{r}\cdot\hat{z}=0\)) and \(\hat{u}=\hat{z}\times\hat{r}\) a
unit vector perpendicular to it in the plane of the screen.
Imagine now going from the observer to some point along a line on
a first screen. It is easiest to work in terms of angles relative
to the observer, so we use \(\rho=p/d\) and \(\varsigma=s/d\) to write this
trajectory as,
(2)\[ d(\hat{z} + \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}).\]
When it hits the line, light can be bent only perpedicular to the
line, by an angle which we will label \(\alpha\) (with positive \(\alpha\) implying
bending closer to the line of sight; equal to \(\hat\alpha\) in Simard &
Pen). Hence, beyond the screen it will travel along
(3)\[ d(\hat{z} + \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}) - (d-d_{1})\alpha_{1}\hat{r}_{1}.\]
If it then hits a line on the second screen, it will again be bent,
and then follow,
(4)\[ d(\hat{z} + \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}) - (d-d_{1})\alpha_{1}\hat{r}_{1} - (d-d_{2})\hat{r}_{2}.\]
In order to specify the full trajectory, we need to make sure that it
actually intersects the line on the second screen, and ends at the
pulsar, i.e.,
(5)\[\begin{split} \begin{eqnarray}
d_{2}(\hat{z} + \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}) - (d_{2}-d_{1})\alpha_{1}\hat{r}_{1}
&=& d_{2}(\hat{z} + \rho_{2}\hat{r}_{2} + \varsigma_{2}\hat{u}_{2}),\\
d_{p}(\hat{z} + \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}) - (d_{p}-d_{1})\alpha_{1}\hat{r}_{1} - (d_{p}-d_{2})\alpha_{2}\hat{r}_{2}
&=& d_{p}\hat{z}.
\end{eqnarray}\end{split}\]
This can be simplified to,
(6)\[\begin{split} \begin{eqnarray}
\rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1} - (1-d_{1}/d_{2})\alpha_{1}\hat{r}_{1}
&=& \rho_{2}\hat{r}_{2} + \varsigma_{2}\hat{u}_{2},\\
\rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1} - (1-d_{1}/d_{p})\alpha_{1}\hat{r}_{1} - (1-d_{2}/d_{p})\hat{r}_{2}) &=& 0.
\end{eqnarray}\end{split}\]
Obviously, this can be extended to multiple screens, and for \(n\)
screens one generally winds up with \(n\) equations for two-dimensional
vectors with \(2n\) unknowns, the \(n\) bending angles \(\alpha_{i}\) and the \(n\)
angular offsets \(\varsigma_{i}\) along the lines.
Direct solutions for one or two screens
If there is only a single screen, the situation is simple: we just
need to make sure the trajectory after the first screen ends at the
pulsar. Following the above simplification, one needs,
(7)\[ \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1} - (1-d_{1}/d_{p})\alpha_{1}\hat{r}_{1} = 0.\]
Thus, taking inner products with \(\hat{r}_{1}\) and \(\hat{u}_{1}\), one finds
\(\alpha_{1}=\rho_{1}d_{p}/(d_{p}-d_{1})\) and \(\varsigma_{1}=0\), i.e., the ray will always be
bent along the point closes to the line of sight.
For two screens, the above still can be solved fairly easily by
substitution, though it becomes less obvious this is still useful.
We start by considering inner products with \(\hat{r}_{1}\)
and \(\hat{u}_{1}\), giving the following four equations (writing
\(s_{12}\equiv1-d_{1}/d_{2}\), etc., for brevity),
\[\begin{split}\begin{eqnarray}
\rho_{1} - s_{12}\alpha_{1} &=& \rho_{2}\hat{r}_{2}\cdot\hat{r}_{1} + \varsigma_{2}\hat{u}_{2}\cdot\hat{r}_{1}
= \rho_{2}\cos \delta - \varsigma_{2}\sin \delta,\\
\varsigma_{1} &=& \rho_{2}\hat{r}_{2}\cdot\hat{u}_{1} + \varsigma_{2}\hat{u}_{2}\cdot\hat{u}_{1}
= \rho_{2} \sin \delta + \varsigma_{2} \cos \delta,\\
\rho_{1} - s_{1p}\alpha_{1} &=& s_{2p}\alpha_{2}\hat{r}_{2}\cdot\hat{r}_{1}
= s_{2p}\alpha_{2}\cos \delta\\
\varsigma_{1} &=& s_{2p}\alpha_{2}\hat{r}_{2}\cdot\hat{u}_{1}
= s_{2p}\alpha_{2}\sin \delta,
\end{eqnarray}\end{split}\]
where in the second equalities we used
\(\hat{r}_{2}\cdot\hat{r}_{1}=\hat{u}_{2}\cdot\hat{u}_{1}=\cos \delta\) and
\(\hat{r}_{2}\cdot\hat{u}_{1}=-\hat{u}_{2}\cdot\hat{r}_{1}=\sin \delta\).
Eliminating \(s_{2p}\alpha_{2}\) between the third and fourth,
\[\rho_{1} - s_{1p}\alpha_{1} = \varsigma_{1}/\tan \delta.\]
Solving for \(\alpha_{1}\) and inserting this in the first,
\[\rho_{1} - (s_{12}/s_{1p})(\rho_{1}-\varsigma_{1}/\tan \delta) = \rho_{2}\cos \delta - \varsigma_{2}\sin \delta.\]
Collecting \(\rho_{1}\) terms and inserting the second,
\[\rho_{1}(1 - s_{12}/s_{1p}) + (s_{12}/s_{1p})(\rho_{2} \sin \delta + \varsigma_{2} \cos \delta)/\tan \delta
= \rho_{2}\cos \delta - \varsigma_{2}\sin \delta.\]
Bringing \(\varsigma_{2}\) and \(\rho_{2}\) terms together, this simplifies to
\[\varsigma_{2}(\sin \delta + (s_{12}/s_{1p}) \cos^{2}\delta/\sin \delta)
= (\rho_{2}\cos \delta - \rho_{1})(1-s_{12}/s_{1p}),\]
and thus
\[\varsigma_{2}(1 - \cos^{2}\delta(1-s_{12}/s_{1p})) = (\rho_{2}\cos \delta -\rho_{1})(1-s_{12}/s_{1p})\sin \delta.\]
The other unknowns then follow from substitution.
General solution
For \(n\) screens, it has to hold for all screens \(i>1\), that,
(8)\[ \rho_{1}\hat{r}_{1} + \varsigma_{1}\hat{u}_{1}
- \sum_{j=1}^{i-1} \alpha_{j}(1-d_{j}/d_{i})\hat{r}_{j}
= \rho_{i}\hat{r}_{i} + \varsigma_{i}\hat{u}_{i}.\]
The same has to hold for the pulsar, and one can incorporate this by
giving it an index \(p\) for which \(p-1=n\) in the summation. Above, we
effectively assumed \(\rho_{p}=\varsigma_{p}=0\), but in the general case one
will have a spatial offset \(\vec{r}_{p}\) corresponding to an angular
offset \(\rho_{p}=r_{p}/d_{p}\) in direction \(\hat{r}_{p}\) (and \(\varsigma_{p}=0\)).
One can also have the telescope be offset (e.g., from the solar system
barycentre), by some \(\vec{r}_{t}\) with norm \(r_{t}\) and direction
\(\hat{r}_{t}\). Essentially the same equation will hold, with summation
starting at \(j=t=0\), \(d_{t}=0\), \(\rho_{ti}=r_{t}/d_{i}\), and \(\alpha_{t}\) and \(\varsigma_{t}\) angles towards
the line of sight and perpendicular to \(\hat{r}_{t}\), respectively.
Bringing unknowns to one side and the knowns to the other, and again
writing \(s_{ji}=1-d_{j}/d_{i}\) (where all \(s_{ti}=1\)), one finds,
(9)\[- (\varsigma_{i}\hat{u}_{i} - \varsigma_{t}\hat{u}_{t})
- \sum_{j=0}^{i-1} \alpha_{j}s_{ji}\hat{r}_{j}
= \rho_{i}\hat{r}_{i} - \rho_{t}\hat{r}_{t}\]
To solve this, we project on given \(x\) and \(y\) directions, i.e., use
\(\hat{r}_{x,y}\) and \(\hat{u}_{x,y}\) (in terms of angles
\(\phi_{i}\) relative to \(\hat{x}\),
one has \(\hat{r}_{i} = \cos \phi_{i} \hat{x} + \sin \phi_{i}\hat{y}\)
and thus \(\hat{u}_{i} = -\sin \phi_{i} \hat{x} + \cos \phi_{i}\hat{y}\)).
Furthermore, for brevity, we write
\(\vec{\theta} \equiv (p_{i}\hat{r}_{i} - p_{t}\hat{r}_{t}) / d_{i}\).
With that, the equations in matrix form are,
(10)\[\begin{split}\left(\begin{matrix}
\hat{u}_{t,x} & -\hat{r}_{t,x} & -\hat{u}_{1,x} & 0 & \ldots & 0 & 0\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & -\hat{u}_{1,y} & 0 & \ldots & 0 & 0\\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{12}\hat{r}_{1,x} & \ldots & 0 & 0\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{12}\hat{r}_{1,y} & \ldots & 0 & 0\\
\vdots & \vdots &\vdots & \vdots & \ddots & \vdots &\vdots \\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{1n}\hat{r}_{1,x} & \ldots & -\hat{u}_{n,x} & 0\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{1n}\hat{r}_{1,y} & \ldots & -\hat{u}_{n,y} & 0\\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{1p}\hat{r}_{1,x} & \ldots & 0 & -s_{np}\hat{r}_{n,x}\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{1p}\hat{r}_{1,y} & \ldots & 0 & -s_{np}\hat{r}_{n,y}\\
\end{matrix}\right)
\left(\begin{matrix}
\varsigma_{t}\\
\alpha_{t}\\
\varsigma_{1}\\
\alpha_{1}\\
\vdots\\
\varsigma_{n-1}\\
\alpha_{n-1}\\
\varsigma_{n}\\
\alpha_{n}
\end{matrix}\right)
= \left(\begin{matrix}
\theta_{1,x}\\
\theta_{1,y}\\
\theta_{2,x}\\
\theta_{2,y}\\
\vdots\\
\theta_{n,x}\\
\theta_{n,y}\\
\theta_{p,x}\\
\theta_{p,y}\\
\end{matrix}\right).\end{split}\]
These can be solved by by inverting the matrix \(A\), i.e.,
(11)\[\begin{split} \left(\begin{matrix}
\varsigma_{t}\\
\alpha_{t}\\
\varsigma_{1}\\
\alpha_{1}\\
\vdots\\
\varsigma_{n-1}\\
\alpha_{n-1}\\
\varsigma_{n}\\
\alpha_{n}
\end{matrix}\right) = A^{-1}
\left(\begin{matrix}
\theta_{1,x}\\
\theta_{1,y}\\
\theta_{2,x}\\
\theta_{2,y}\\
\vdots\\
\theta_{n,x}\\
\theta_{n,y}\\
\theta_{p,x}\\
\theta_{p,y}\\
\end{matrix}\right).\end{split}\]
Velocities
In principle, the telescope, screens and pulsar will all have velocities. One
sees that the entries in the matrix involve only directions of the telescope and
screens and ratios of distances, which will change with time much more slowly
than everything else. Hence, one can solve for the time derivatives of the
parameters by applying the matrix inverse to the time derivatives of the entries
of the right-hand side vector, which are simply the \(x\) and \(y\)
components of the proper motions of the screens and the pulsar relative to the
telescope, i.e.,
(12)\[\begin{split}\left(\begin{matrix}
\dot{\varsigma}_{t}\\
\dot{\alpha}_{t}\\
\dot{\varsigma}_{1}\\
\dot{\alpha}_{1}\\
\vdots\\
\dot{\varsigma}_{n-1}\\
\dot{\alpha}_{n-1}\\
\dot{\varsigma}_{n}\\
\dot{\alpha}_{n}
\end{matrix}\right)
= A^{-1}
\left(\begin{matrix}
\mu_{1,x}\\
\mu_{1,y}\\
\mu_{2,x}\\
\mu_{2,y}\\
\vdots\\
\mu_{n,x}\\
\mu_{n,y}\\
\mu_{p,x}\\
\mu_{p,y}\\
\end{matrix}\right).\end{split}\]
Trajectories and time delays
The above \(\theta\) and \(\mu\) reflect the positions and proper
motions of the structures, not of trajectories taken by different rays.
Those are given by,
(13)\[\begin{split} \begin{eqnarray}
\vec{\theta}^{r}_{i} &=& \theta^{s}_{i}\hat{r}_{i} + \varsigma_{i} \hat{u}_{i},\\
\vec{\mu}^{r}_{i} &=& \mu^{s}_{i}\hat{r}_{i} + \dot{\varsigma}_{i} \hat{u}_{i}.
\end{eqnarray}\end{split}\]
The total time delay \(\tau\) and its derivative \(\dot{\tau}\) are given
by,
(14)\[\begin{split} \begin{eqnarray}
\tau &=& \sum_{i=0}^{n} \frac{d_{i+1}-d_{i}}{2c}
\left|\vec{\theta}^{r}_{i+1}-\vec{\theta}^{r}_{i}\right|^{2},\\
\dot{\tau} &=&
\sum_{i=0}^{n} \frac{d_{i+1}-d_{i}}{c}
\left(\vec{\theta}^{r}_{i+1}-\vec{\theta}^{r}_{i}\right)
\cdot\left(\vec{\mu}^{r}_{i+1}-\vec{\mu}^{r}_{i}\right).
\end{eqnarray}\end{split}\]
Bending angle dependent positions
In the derivation above, the rays have to intersect the one dimensional lines on
each screen. Physically, this corresponds to assuming the lensing structures
have widths much smaller than their separations. In reality, where a ray
crosses through a lens will depend on the bending angle: very small angles are
produced right on top of the lens and far away from it, while the largest
bending angles occur where the electron column density gradient is steepest.
While the general bending-angle dependence of the crossing location depends on
the precise lens shape, and thus cannot be easily included, it turns out to be
nearly trivial to include a linear expansion, where the angular location of the
crossing point is given by,
(15)\[p(\alpha) = p(0) + \alpha \frac{dp}{d\alpha}
\quad\Rightarrow\quad
\rho(\alpha) = \rho(0) + \frac{\alpha}{d}\frac{dp}{d\alpha}
= \rho(0) + \alpha \rho^{\prime},\]
where in the last equality we implicitly defined \(\rho^{\prime}\equiv
(1/d)(dp/d\alpha)\).
The case described so far then corresponds to the zeroth-order approximation,
with \(dp/d\alpha=0\).
Looking at the general solution (Eq. (9)), one sees
that the only difference is that the \(\rho_{i}\hat{r}_{i}\) term on the
right now no longer is constant.
Taking the \(\alpha\)-dependent part to the left, one finds,
(16)\[ - (\varsigma_{i}\hat{u}_{i} - \varsigma_{t}\hat{u}_{t})
- \alpha_{i} \rho_{i}^{\prime} \hat{r}_{i}
- \sum_{j=0}^{i-1} \alpha_{j}s_{ji}\hat{r}_{j}
= \rho_{i}(0)\hat{r}_{i} - \rho_{t}\hat{r}_{t},\]
In matrix form, one then has,
(17)\[\begin{split}\left(\begin{matrix}
\hat{u}_{t,x} & -\hat{r}_{t,x} & -\hat{u}_{1,x} & -\rho_{1}^{\prime}\hat{r}_{1,x} & \ldots & 0 & 0\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & -\hat{u}_{1,y} & -\rho_{1}^{\prime}\hat{r}_{1,y} & \ldots & 0 & 0\\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{12}\hat{r}_{1,x} & \ldots & 0 & 0\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{12}\hat{r}_{1,y} & \ldots & 0 & 0\\
\vdots & \vdots &\vdots & \vdots & \ddots & \vdots &\vdots \\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{1n}\hat{r}_{1,x} & \ldots & -\hat{u}_{n,x} & -\rho_{n}^{\prime} \hat{r}_{n,x}\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{1n}\hat{r}_{1,y} & \ldots & -\hat{u}_{n,y} & -\rho_{n}^{\prime} \hat{r}_{n,y}\\
\hat{u}_{t,x} & -\hat{r}_{t,x} & 0 & -s_{1p}\hat{r}_{1,x} & \ldots & 0 & -s_{np}\hat{r}_{n,x}\\
\hat{u}_{t,y} & -\hat{r}_{t,y} & 0 & -s_{1p}\hat{r}_{1,y} & \ldots & 0 & -s_{np}\hat{r}_{n,y}\\
\end{matrix}\right)
\left(\begin{matrix}
\varsigma_{t}\\
\alpha_{t}\\
\varsigma_{1}\\
\alpha_{1}\\
\vdots\\
\varsigma_{n-1}\\
\alpha_{n-1}\\
\varsigma_{n}\\
\alpha_{n}
\end{matrix}\right)
= \left(\begin{matrix}
\theta_{1,x}\\
\theta_{1,y}\\
\theta_{2,x}\\
\theta_{2,y}\\
\vdots\\
\theta_{n,x}\\
\theta_{n,y}\\
\theta_{p,x}\\
\theta_{p,y}\\
\end{matrix}\right).\end{split}\]
Like before, one can use the same matrix inverse to calculate velocities.
The trajectories relative to the telescope now are given by,
(18)\[\begin{split} \begin{eqnarray}
\theta^{r}_{i}
= \left(\theta^{s}_{i} + \alpha_{i}\rho_{i}^{\prime}\right)\hat{r}_{i}
+ \varsigma_{i} \hat{u}_{i},\\
\mu^{r}_{i}
= \left(\mu^{s}_{i} + \dot{\alpha}_{i}\rho_{i}^{\prime}\right)\hat{r}_{i}
+ \dot{\varsigma}_{i} \hat{u}_{i}.
\end{eqnarray}\end{split}\]