Simultaneous GPs on photometry and spectroscopy

Simultaneous GPs on photometry and spectroscopy#

Imagine we have both photometric and spectroscopic time series of a rather active star. Photometry will provide essential information on the rotational period of the star and the decay time scale of activity regions, while spectroscopic activity indexes are essential to disentangle planetary and stellar signal in radial velocities.

In this example, we will analyze the planetary system around K2-141 Malavolta et al. 2018. The spectroscopic dataset is taken from Bonomo et al. 2023 and includes Radial velocity data, Bisector Inverse Span, and S-index, all gathered with HARPS-N.

Photometric data comes from K2 campaign 12, corrected fro systematics as described in Vanderburg and Johnson 2014 and available at the MAST K2SFF website. Data from TESS sector 42 and sector 70 are available, but we are not going to use them in this example.

lightcurve of K2-141

Fig. 4 Photometric data for K2-141. The stellar activity modulation is clearly identifiable.#

We want to model stellar activity in photometry simultaneously with the spectroscopic time series.

  • For photometric data, we use the [exponential-sine periodic kernel] (ESP) implemented in ‘S+LEAF’. Each dataset will have its independent covariance matrix.

  • For spectroscopic data, we use the multidimensional quasi-periodic kernel through tinygp.

For the photometric data, we also include the transit model for the two planets, and a normalization factor, independent for each dataset (if you have more than one light curve). Check the Light curve fit section in Quickstart for additional information.

For the spectroscopic data, we only have to include the radial_velocities model due to the planets in the radial velocity dataset.

The input section takes this form:

 1inputs:
 2  LC_K2_c12:
 3    file: datasets/K2-141_K2_KSFF_PyORBIT.dat
 4    kind: Phot
 5    models:
 6      - spleaf_esp
 7      - lc_model
 8      - normalization_factor
 9  RVdata:
10    file: datasets/K2-141_RV_PyORBIT.dat
11    kind: RV
12    models:
13      - radial_velocities
14      - gp_multidimensional
15  BISdata:
16    file: datasets/K2-141_BIS_PyORBIT.dat
17    kind: BIS
18    models:
19      - gp_multidimensional
20  Sdata:
21    file: datasets/K2-141_Sindex_PyORBIT.dat
22    kind: S_index
23    models:
24      - gp_multidimensional

In the common section, we specify:

  • the two transiting planets: planets:b and planet:c. Note that we do not specify any prior on the period or transit time, as they will be fitted directly in photometric data;

  • activity: the rotational period of the star Prot, the decay timescale of activity regions Pdec, and the coherence scale Oamp, will be shared among the trained ESP and the multidimensional QP.

  • stellar parameters: in addition to he priors on mass, radius, and density, we also include priors on the limb darkening coefficients for the two instruments involved in our analysis;

25common:
26  planets:
27    b:
28      orbit: circular
29      use_time_inferior_conjunction: True
30      boundaries:
31        P: [0.2750, 0.2850]
32        K: [0.001, 20.0]
33        Tc: [57744.00, 57744.10]
34      spaces:
35        P: Linear
36        K: Linear
37    c:
38      orbit: keplerian
39      parametrization: Eastman2013
40      use_time_inferior_conjunction: True
41      boundaries:
42        P: [7.70, 7.80]
43        K: [0.001, 20.0]
44        Tc: [58371.00, 58371.10]
45        e: [0.00, 0.70]
46      priors:
47        e: ['Gaussian', 0.00, 0.098]
48      spaces:
49        P: Linear
50        K: Linear
51  activity:
52    boundaries:
53      Prot: [10.0, 20.0]
54      Pdec: [20.0, 1000.0]
55      Oamp: [0.001, 1.0]
56  star:
57    star_parameters:
58      priors:
59        mass: ['Gaussian', 0.708, 0.028]
60        radius: ['Gaussian', 0.681, 0.018]
61        density: ['Gaussian', 2.65, 0.08]
62    limb_darkening:
63      type: ld_quadratic
64      priors:
65        ld_c1: ['Gaussian', 0.68, 0.10]
66        ld_c2: ['Gaussian', 0.05, 0.10]

In the models section we include two models for Gaussian Process regression: the spleaf_esp will perform GP regression on each light curve independently (if you have more than one light curve), i.e., each light curve will have its own covariance matrix and its own amplitude of the covariance, similarly when performing the trained GP regression.

The gp_multidimensional will perform multidimensional GP regression over the datasets that are employing these model, namely for this specific example, the RV dataset, the BIs dataset, and the S index.

You may notice that both spleaf_esp and gp_multidimensional are recalling the same common model common:activity. The GP hyperparameters will be shared, in this sense we may say that photometry is training the hyperparameters of the multidimensional GP.

 67models:
 68  radial_velocities:
 69    planets:
 70      - b
 71      - c
 72  lc_model:
 73    kind: pytransit_transit
 74    limb_darkening: limb_darkening
 75    planets:
 76      - b
 77      - c
 78  normalization_factor:
 79    kind: local_normalization_factor
 80    boundaries:
 81      n_factor: [0.9, 1.1]
 82  spleaf_esp:
 83    model: spleaf_esp
 84    common: activity
 85    n_harmonics: 4
 86    hyperparameters_condition: True
 87    rotation_decay_condition: True
 88    LC_K2_c12:
 89      boundaries:
 90        Hamp: [0.000, 1.00]
 91  gp_multidimensional:
 92    model: tinygp_multidimensional_quasiperiodic
 93    common: activity
 94    hyperparameters_condition: True
 95    rotation_decay_condition: True
 96    RVdata:
 97      boundaries:
 98        rot_amp: [-100.0, 100.0] #at least one must be positive definite
 99        con_amp: [0.0, 100.0]
100      derivative: True
101    BISdata:
102      boundaries:
103        rot_amp: [-100.0, 100.0]
104        con_amp: [-100.0, 100.0]
105      derivative: True
106    Sdata:
107      boundaries:
108        con_amp: [-1.0, 1.0]
109      derivative: False

The rest of the configuration file looks as usual.

129parameters:
130  Tref: 59200.00
131  safe_reload: True
132  use_tex: False
133  low_ram_plot: True
134  plot_split_threshold: 1000
135  cpu_threads: 32
136solver:
137  pyde:
138    ngen: 50000
139    npop_mult: 4
140  emcee:
141    npop_mult: 4
142    nsteps: 50000
143    nburn: 20000
144    nsave: 25000
145    thin: 100
146  recenter_bounds: True

Sharing only some of the hyperparameters#

Improved interface starting from version 10.9

The following example applies to analysis performed with PyORBIT version >= 10.9 Previous version salready implemented the following analysis, with a less user-friendly interface.

The previous configuration assumes that you want to share the rotational period of the star Prot, the decay time scale of active regions Pdec, and the coherence scale Oamp among all the datasets. There are reasons to believe that the Pdec may change with time, and Oamp may change from datasets to dataset, depending on the specific aspects of stellar activity captured by each observational method. To share only some of the hyperparameters involved in the GP regression, we associated the photometric and the spectroscopic models to two separate common activity models. In this way, the hyperparameters will be independent. Since it would not make sense to have all the hyperparameters being independent (as it wold be equivalent to run two independent fits), we specify that the rotational period of the star Prot will be extracted from the stellar properties rather then from the activity model by setting the use_stellar_rotation_period flag to True:

  activity_photometry:
    model: activity
    use_stellar_rotation_period: True
    boundaries:
      #Prot: [10.0, 20.0]
      Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]
  activity_spectroscopy:
    model: activity
    use_stellar_rotation_period: True
    boundaries:
      #Prot: [10.0, 20.0]
      Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]

In this way, there will be independent Pdec and Oamp depending on which common activity model each model will recall.

Note how we have to specify two different names for the models, activity_photometry and activity_spectroscopy, while declaring that they are both model:activity. Since the rotational period hyperparameter is not taken from the activity model anymore, there is no need to specify its boundaries in the activity model. Instead , you can specifiy boundaries, prior, and space foir the rotation_period in the star_parameters section.

  star:
    star_parameters:
      boundaries:
        rotation_period: [10.0, 20.0]
      priors:
        mass: ['Gaussian', 0.708, 0.028]
        radius: ['Gaussian', 0.681, 0.018]
        density: ['Gaussian', 2.65, 0.08]
    limb_darkening:
      type: ld_quadratic
      priors:
        ld_c1: ['Gaussian', 0.68, 0.10]
        ld_c2: ['Gaussian', 0.05, 0.10]
   #use_stellar_activity_decay
   #activity_decay

Finally, remember to associate the corresponding common activity to the model used for each dataset, as in the example below (with some omessions to highlight the differences with respect to the previous code block).

models:
  spleaf_esp:
    model: spleaf_esp
    common: activity_photometry 
    n_harmonics: 4
    hyperparameters_condition: True
    ...
  gp_multidimensional:
    model: tinygp_multidimensional_quasiperiodic
    common: activity_spectroscopy
    hyperparameters_condition: True
    ...

The rotation_period parameter will be taken automatically from the stellar porperties. If there is only one stellar model, the association will be automatic without requiring the explicit declaration in the models:*model_name*:common section.

If you want to share also the decay time scale Pdec, you can follow the same procedure as the rotational period Prot, noting that the flag to be activated is use_stellar_activity_decay rather than use_stellar_rotation_period. In the example below, both the Prot and Pdec are shared among the models.

  activity_photometry:
    model: activity
    use_stellar_rotation_period: True
    use_stellar_activity_decay: True
    boundaries:
      #Prot: [10.0, 20.0]
      #Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]
  activity_spectroscopy:
    model: activity
    use_stellar_rotation_period: True
    use_stellar_activity_decay: True
    boundaries:
      #Prot: [10.0, 20.0]
      #Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]

In the stellar properties, the decay time scale takes the name of activity_decay:

  star:
    star_parameters:
      boundaries:
        rotation_period: [10.0, 20.0]
        activity_decay: [20.0, 1000.0]
      priors:
        mass: ['Gaussian', 0.708, 0.028]
        ...

The models ection would be the same as in the case of the shared rotational period.

Finally, it is possible to create mixed cases where a given hyperparameters is shared among some datasets but not with others. In the following (fictional) examples, the Prot is shared among all the parameters, while `Pdec’ is shared only among the photometric models.

  activity_ground_photometry:
    model: activity
    use_stellar_rotation_period: True
    use_stellar_activity_decay: True
    boundaries:
      #Prot: [10.0, 20.0]
      #Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]
  activity_space_photometry:
    model: activity
    use_stellar_rotation_period: True
    use_stellar_activity_decay: True
    boundaries:
      #Prot: [10.0, 20.0]
      #Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]
activity_spectroscopy:
    model: activity
    use_stellar_rotation_period: True
    boundaries:
      #Prot: [10.0, 20.0]
      Pdec: [20.0, 1000.0]
      Oamp: [0.001, 1.0]

Note that omitting the declaration for use_stellar_activity_decay is equivalent to declare use_stellar_activity_decay: False.

Tip

For now, there is no flag to share the coherence scale Oamp, as this value is the last one you want to share among models. In other words, if you inclided to share Oamp, then likely you wantto share Prot and Pdec as well, following back to the case highlighted at the beginning of this page.