.

Rossiter-McLaughlin Effect¶

Setup¶

%matplotlib inline


As always, let’s do imports and initialize a logger and a new bundle. See Building a System for more details.

import phoebe
from phoebe import u # units
import numpy as np
import matplotlib.pyplot as plt

logger = phoebe.logger()

b = phoebe.default_binary()


Let’s make a significant mass ratio and radius ratio...

b['q'] = 0.7

b['requiv@primary'] = 1.0

b['requiv@secondary'] = 0.5

b['teff@secondary@component'] = 5000


Make sure the primary star is spinning quickly...

b['syncpar@primary@component'] = 2


Now we’ll add radial velocity and mesh datasets. We’ll add two identical datasets for RVs so that we can have one computed dynamically and the other computed numerically (these options will need to be set later).

b.add_dataset('rv', times=np.linspace(0,2,201), dataset='dynamicalrvs')

<ParameterSet: 15 parameters | contexts: compute, dataset>

b.add_dataset('rv', times=np.linspace(0,2,201), dataset='numericalrvs')

<ParameterSet: 15 parameters | contexts: compute, dataset>


Storing the mesh at every timestep is overkill, and will be both computationally and memory intensive. So let’s just sample at the times we care about.

times = b.get_value('times@primary@numericalrvs@dataset')
times = times[times<0.1]
print times

[0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09]

b.add_dataset('mesh', dataset='mesh01', times=times, columns=['vws', 'rvs*'])

<ParameterSet: 4 parameters | contexts: compute, dataset>


Running Compute¶

Now let’s set the rv_method so that one dataset uses the dynamical method and the other uses the flux-weighted (numerical) method. Note that here we have to use set_value_all or loop over the components, as technically there are parameters for each component-dataset pair.

b.set_value_all('rv_method@dynamicalrvs@compute', 'dynamical')
b.set_value_all('rv_method@numericalrvs@compute', 'flux-weighted')


Let’s check to make sure that rv_method is set as we’d expect.

print b['rv_method']

ParameterSet: 4 parameters
rv_method@primary@dynamical...: dynamical
rv_method@secondary@dynamic...: dynamical
rv_method@primary@numerical...: flux-weighted
rv_method@secondary@numeric...: flux-weighted

b.run_compute(irrad_method='none')

<ParameterSet: 110 parameters | kinds: rv, mesh>


Plotting¶

Now let’s plot the radial velocities.

First we’ll plot the dynamical RVs. Note that dynamical RVs show the true radial velocity of the center of mass of each star, and so we do not see the Rossiter McLaughlin effect.

afig, mplfig = b['dynamicalrvs@model'].plot(c={'primary': 'b', 'secondary': 'r'}, show=True)


But the numerical method integrates over the visible surface elements, giving us what we’d observe if deriving RVs from observed spectra of the binary. Here we do see the Rossiter McLaughlin effect. You’ll also notice that RVs are not available for the secondary star when its completely occulted (they’re nans in the array).

afig, mplfig = b['numericalrvs@model'].plot(c={'primary': 'b', 'secondary': 'r'}, show=True)


To visualize what is happening, we can plot the radial velocities of each surface element in the mesh at one of these times.

Here just plot on the mesh@model parameterset - the mesh will automatically get coordinates from mesh01 and then we point to rvs@numericalrvs for the facecolors.

afig, mplfig = b['mesh@model'].plot(time=0.03, fc='rvs@numericalrvs', ec="None", show=True)


Here you can see that the secondary star is blocking part of the “red” RVs of the primary star.

This is essentially the same as plotting the negative z-component of the velocities (for convention - our system is in a right handed system with +z towards the viewer, but RV convention has negative RVs for blue shifts).

We could also plot the RV per triangle by plotting ‘vws’. Note that this is actually defaulting to an inverted colormap to show you the same colorscheme (‘RdBu_r’ vs ‘RdBu’).

afig, mplfig = b['mesh01@model'].plot(time=0.09, fc='vws', ec="None", show=True)

.
Last update: 10/29/2018 9:20 a.m. (CET)