1D RM-Synthesis¶

[1]:

from __future__ import annotations

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.visualization import quantity_support
from rm_lite.tools_1d import rmsynth

plt.rcParams["figure.dpi"] = 150

_ = quantity_support()
rng = np.random.default_rng(42)

First generate some synthetic data

[2]:

from rm_lite.utils.simulate import faraday_simple_spectrum, faraday_slab_spectrum
from rm_lite.utils.synthesis import freq_to_lambda2

Here we’ll simulate RACS-all frequency coverage

[3]:

bw_low = 288
bw_mid = 144
bw_high = 288
low = np.linspace(943.5 - bw_low / 2, 943.5 + bw_low / 2, 36) * u.MHz
mid = np.linspace(1367.5 - bw_mid / 2, 1367.5 + bw_mid / 2, 9) * u.MHz
high = np.linspace(1655.5 - bw_high / 2, 1655.5 + bw_high / 2, 9) * u.MHz
freqs = np.concatenate([low, mid, high])
freq_hz = freqs.to(u.Hz).value

Now we make a Faraday simple spectrum with a single RM component. We will use the following parameters:

[4]:

delta_rm_radm2 = 30
rm_radm2 = 100
frac_pol = 0.5
psi0_deg = 10
complex_data_noiseless = faraday_simple_spectrum(
    freq_to_lambda2(freq_hz),
    frac_pol=frac_pol,
    psi0_deg=psi0_deg,
    rm_radm2=rm_radm2,
)

[5]:

fig, ax = plt.subplots()
ax.plot(
    freq_hz, np.real(complex_data_noiseless), ".", label="Stokes Q", color="tab:red"
)
ax.plot(
    freq_hz, np.imag(complex_data_noiseless), ".", label="Stokes U", color="tab:blue"
)
ax.legend()
ax.set(
    xlabel=rf"$\nu$ / {u.Hz:latex_inline}",
    ylabel="Flux density",
    title="Stokes Q and U",
)

[5]:

[Text(0.5, 0, '$\\nu$ / $\\mathrm{Hz}$'),
 Text(0, 0.5, 'Flux density'),
 Text(0.5, 1.0, 'Stokes Q and U')]

Now we can run RM-synthesis by calling rmsynth.run_rmsynth

[6]:

help(rmsynth.run_rmsynth)

Help on function run_rmsynth in module rm_lite.tools_1d.rmsynth:

run_rmsynth(freq_arr_hz: 'NDArray[np.float64]', complex_pol_arr: 'NDArray[np.complex128]', complex_pol_error: 'NDArray[np.complex128]', stokes_i_arr: 'NDArray[np.float64] | None' = None, stokes_i_error_arr: 'NDArray[np.float64] | None' = None, stokes_i_model_arr: 'NDArray[np.float64] | None' = None, stokes_i_model_error: 'NDArray[np.float64] | None' = None, phi_max_radm2: 'float | None' = None, d_phi_radm2: 'float | None' = None, n_samples: 'float | None' = 10.0, weight_type: 'WeightType' = 'variance', robust: 'float | None' = None, do_fit_rmsf: 'bool' = False, do_fit_rmsf_real: 'bool' = False, fit_function: "Literal['log', 'linear']" = 'log', fit_order: 'int' = 2, ignore_stokes_i: 'bool' = False, moment_threshold_snr: 'float' = 5.0) -> 'RMSynth1DResults'
    Run RM-synthesis on 1D data

    Args:
        freq_arr_hz (NDArray[np.float64]): Frequencies in Hz
        complex_pol_arr (NDArray[np.complex128]): Complex polarisation values (Q + iU)
        complex_pol_error (NDArray[np.float64]): Complex polarisation errors (dQ + idU)
        stokes_i_arr (NDArray[np.float64] | None, optional): Total itensity values. Defaults to None.
        stokes_i_error_arr (NDArray[np.float64] | None, optional): Total intensity errors. Defaults to None.
        stokes_i_model_arr (NDArray[np.float64] | None, optional): Total intensity model array. Defaults to None.
        stokes_i_model_error (NDArray[np.float64] | None, optional): Total intensity model error. Defaults to None.
        phi_max_radm2 (float | None, optional): Maximum Faraday depth. Defaults to None.
        d_phi_radm2 (float | None, optional): Spacing in Faraday depth. Defaults to None.
        n_samples (float | None, optional): Number of samples across the RMSF. Defaults to 10.0.
        weight_type (WeightType, optional): Weighting: 'variance' (1/sigma^2), 'uniform' (equal per channel), 'uniform_lsq' (equal per lambda^2 interval, narrows the RMSF), 'briggs' (robust). Defaults to "variance".
        robust (float | None, optional): Briggs robust parameter, required for weight_type='briggs'. Defaults to None.
        do_fit_rmsf (bool, optional): Fit the RMSF main lobe. Defaults to False.
        do_fit_rmsf_real (bool, optional): Fit only the real part of the RMSF. Defaults to False.
        fit_function ("log" | "linear", optional): RMSF fit function. Defaults to "log".
        fit_order (int, optional): Polynomial fit order. Defaults to 2. Negative values will iterate until the fit is good.
        moment_threshold_snr (float, optional): SNR cut (times the theoretical FDF noise) applied to FDF amplitudes before computing the Faraday moments. Defaults to 5.0.

    Returns:
        RMSynth1DResults:
            fdf_parameters (pl.DataFrame): FDF parameters
            fdf_arrs (pl.DataFrame): RMSynth arrays
            rmsf_arrs (pl.DataFrame): RMSF arrays

[7]:

fdf_parameters, fdf_arrs, rmsf_arrs, stokes_i_arrs = rmsynth.run_rmsynth(
    freq_arr_hz=freq_hz,
    complex_pol_arr=complex_data_noiseless,
    complex_pol_error=np.zeros_like(complex_data_noiseless),
    do_fit_rmsf=True,
    n_samples=100,
)

WARNING rmsynth.run_rmsynth: Stokes I array/errors or model not provided. No fractional polarization will be calculated.
INFO synthesis.rmsynth_nufft: Running RM-synthesis using the NUFFTs over 2001 Faraday depth channels.
INFO synthesis.rmsynth_nufft: NUFFT complete in 0.0038 seconds.
INFO synthesis.get_rmsf_nufft: Fitting main lobe in each RMSF spectrum.
INFO rmsynth._run_rmsynth: RM-synthesis completed in 14.16ms.

The output values are Polars dataframes that can be inspected easily

[8]:

fdf_parameters

[8]:

shape: (1, 36)

fdf_error_mad	peak_pi_fit	peak_pi_error	peak_pi_fit_debias	peak_pi_fit_snr	peak_pi_fit_index	peak_rm_fit	peak_rm_fit_error	peak_q_fit	peak_u_fit	peak_pa_fit_deg	peak_pa_fit_deg_error	peak_pa0_fit_deg	peak_pa0_fit_deg_error	fit_function	lam_sq_0_m2	ref_freq_hz	fwhm_rmsf_radm2	phi_max_scale_radm2	fdf_error_noise	fdf_q_noise	fdf_u_noise	min_freq_hz	max_freq_hz	n_channels	median_d_freq_hz	frac_pol	frac_pol_error	sigma_add	sigma_add_minus	sigma_add_plus	mom0	mom0_debias	mom1_radm2	mom2_radm2	moment_threshold_snr
f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	str	f64	f64	f64	f64	f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64
NaN	0.500279	0.0	0.500279	inf	1296	99.999904	0.0	-0.202569	-0.457124	123.050018	0.0	10.000452	0.0	"log"	0.082563	1.0433e9	34.522155	113.191071	0.0	0.0	0.0	7.995e8	1.7995e9	54	8.2286e6	NaN	NaN	NaN	NaN	NaN	1.346747	1.346747	69.386287	138.885199	5.0

[9]:

fdf_arrs

[9]:

shape: (2_001, 2)

phi_arr_radm2	fdf_dirty_complex_arr
f64	object
-336.725921	(-0.012941952987807103+0.005575905271668905j)
-336.389195	(-0.013584447466549732+0.005189884663585659j)
-336.052469	(-0.014213132463713843+0.004794461009568579j)
-335.715743	(-0.01482712771337132+0.004389828120960796j)
-335.379017	(-0.015425552941836488+0.003976168117102392j)
…	…
335.379017	(-0.0009008057626928257+0.008053597616970908j)
335.715743	(-0.0009576261370891269+0.007885479952290737j)
336.052469	(-0.0009822552939409935+0.007717971936861509j)
336.389195	(-0.0009751402168175126+0.007551392504192875j)
336.725921	(-0.0009367721104613681+0.007386049740302505j)

[10]:

rmsf_arrs

[10]:

shape: (4_003, 2)

phi2_arr_radm2	rmsf_complex_arr
f64	object
-673.788568	(-0.0329575974415113-0.0360831117177476j)
-673.451842	(-0.032467101832120024-0.03631435677037289j)
-673.115116	(-0.03194042093852627-0.036538602915399204j)
-672.77839	(-0.03137760921834162-0.03675582151570247j)
-672.441664	(-0.030778770230473887-0.03696597750309249j)
…	…
672.441664	(-0.030778770230891164+0.036965977502953294j)
672.77839	(-0.031377609218734684+0.036755821515558944j)
673.115116	(-0.031940420938894855+0.03653860291525102j)
673.451842	(-0.03246710183246473+0.03631435677022105j)
673.788568	(-0.03295760578020578+0.03608308730681604j)

Since we provided no Stokes \(I\) data, the stokes I model will just be unity with 0 error. The flag_arr array tells us which channels were not used in RM-synthesis or model fitting

[11]:

stokes_i_arrs

[11]:

shape: (54, 7)

freq_arr_hz	lambda_sq_arr_m2	stokes_i_model_arr	stokes_i_model_error	flag_arr	complex_pol_arr	complex_pol_error
f64	f64	f64	f64	bool	object	object
7.995e8	0.140606	null	null	true	(-0.49042944532403415-0.09735994638022451j)	0j
8.0773e8	0.137756	null	null	true	(-0.4654241884816873+0.18270283187778674j)	0j
8.1596e8	0.134992	null	null	true	(-0.3001379571186853+0.3998964949791661j)	0j
8.2419e8	0.13231	null	null	true	(-0.053617325136419355+0.497116870006657j)	0j
8.3241e8	0.129707	null	null	true	(0.20074123829142537+0.45793335240974226j)	0j
…	…	…	…	…	…	…
1.6555e9	0.032793	null	null	true	(0.4056251419940925+0.29234952399871j)	0j
1.6915e9	0.031412	null	null	true	(0.46997598588767076+0.17065336998990638j)	0j
1.7275e9	0.030117	null	null	true	(0.49801253407054746+0.04453668048509923j)	0j
1.7635e9	0.0289	null	null	true	(0.49406601726750204-0.07680345422849415j)	0j
1.7995e9	0.027755	null	null	true	(0.463743096038396-0.1869287053309979j)	0j

We can also easily visualise the data

[12]:

phi_arr_radm2 = fdf_arrs["phi_arr_radm2"].to_numpy()
fdf_dirty_arr = fdf_arrs["fdf_dirty_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
x1, x2, y1, y2 = 95, 105, 0.45, 0.55  # subregion of the original image
axins = ax.inset_axes(
    (0.9, 0.6, 0.4, 0.4), xlim=(x1, x2), ylim=(y1, y2), xticklabels=[], yticklabels=[]
)
for _ax in [ax, axins]:
    _ax.plot(
        phi_arr_radm2,
        fdf_dirty_arr.real,
        color="tab:red",
        label="Stokes Q",
    )
    _ax.plot(
        phi_arr_radm2,
        fdf_dirty_arr.imag,
        color="tab:blue",
        label="Stokes U",
    )
    _ax.plot(
        phi_arr_radm2,
        np.abs(fdf_dirty_arr),
        color="k",
        label="Polarized intensity",
    )

    _ax.errorbar(
        fdf_parameters["peak_rm_fit"],
        fdf_parameters["peak_pi_fit"],
        xerr=fdf_parameters["peak_rm_fit_error"],
        yerr=fdf_parameters["peak_pi_error"],
        fmt="o",
        lw=1,
        color="red",
        mfc="none",
        label="Fitted peak",
    )

ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="Flux density",
    title="Dirty FDF",
    # xlim=[50, 150],
)
ax.indicate_inset_zoom(axins, edgecolor="black")
ax.legend()

[12]:

<matplotlib.legend.Legend at 0x7d13a1d29460>

[13]:

phi2_arr_radm2 = rmsf_arrs["phi2_arr_radm2"].to_numpy()
rmsf_arr = rmsf_arrs["rmsf_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi2_arr_radm2,
    np.abs(rmsf_arr),
    color="k",
    label="Polarized intensity",
)
ax.legend()
ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="RMSF",
    title="RMSF",
)

[13]:

[Text(0.5, 0, '$\\phi$ / $\\mathrm{rad\\,m^{-2}}$'),
 Text(0, 0.5, 'RMSF'),
 Text(0.5, 1.0, 'RMSF')]

Now lets do a more complex example. We’ll add noise and a Stokes \(I\) spectrum

[14]:

from rm_lite.utils.fitting import power_law

[15]:

delta_rm_radm2 = 30
rm_radm2 = 100
frac_pol = 0.5
psi0_deg = 10
complex_data_noiseless = faraday_slab_spectrum(
    freq_to_lambda2(freq_hz),
    frac_pol=frac_pol,
    psi0_deg=psi0_deg,
    rm_radm2=rm_radm2,
    delta_rm_radm2=delta_rm_radm2,
)


stokes_i_flux = 1.0
spectral_index = -0.7
rms_noise = 0.1


stokes_i_model = power_law(order=1)
stokes_i_noiseless = stokes_i_model(
    freq_hz / (np.mean(freq_hz)), stokes_i_flux, spectral_index
)
stokes_i_noise = rng.normal(0, rms_noise, size=freq_hz.size)
stokes_i_noisy = stokes_i_noiseless + stokes_i_noise


stokes_q_noise = rng.normal(0, rms_noise, size=freq_hz.size)
stokes_u_noise = rng.normal(0, rms_noise, size=freq_hz.size)
complex_noise = stokes_q_noise + 1j * stokes_u_noise

complex_flux = complex_data_noiseless * stokes_i_noiseless
complex_data_noisy = complex_data_noiseless + complex_noise

Now we enable Stokes \(I\) model fitting through providing the data, and enabling fit_order. If fit_order<0 an iterative fit will be performed.

[16]:

fdf_parameters, fdf_arrs, rmsf_arrs, stokes_i_arrs = rmsynth.run_rmsynth(
    freq_arr_hz=freq_hz,
    complex_pol_arr=complex_data_noisy,
    complex_pol_error=np.ones_like(complex_data_noiseless)
    * (rms_noise + rms_noise * 1j),
    stokes_i_arr=stokes_i_noisy,
    stokes_i_error_arr=np.ones_like(stokes_i_noisy) * rms_noise,
    do_fit_rmsf=True,
    n_samples=100,
    fit_order=-3,
)

INFO synthesis.create_fractional_spectra: Fitting Stokes I model to calculate fractional spectra.
INFO fitting.dynamic_fit: Iteratively fitting Stokes I model of type log with max order 3.
INFO fitting.static_fit: Fitting Stokes I model of type log with order 0.
INFO fitting.static_fit: Fit results: ['1.04 +/- 0.0136']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 1.
INFO fitting.static_fit: Fit results: ['1.07 +/- 0.0138', '-0.671 +/- 0.0608']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 2.
INFO fitting.static_fit: Fit results: ['1.08 +/- 0.0205', '-0.62 +/- 0.0835', '-0.593 +/- 0.672']
INFO fitting.static_fit: Fitting Stokes I model of type log with order 3.
INFO fitting.static_fit: Fit results: ['1.08 +/- 0.0212', '-0.62 +/- 0.127', '-0.589 +/- 1.1', '-0.0254 +/- 6.9']
INFO fitting.dynamic_fit: Fit results for orders [0 1 2 3]:
INFO fitting.dynamic_fit: Best fit found with 2 parameters.
INFO synthesis.rmsynth_nufft: Running RM-synthesis using the NUFFTs over 2001 Faraday depth channels.
INFO synthesis.rmsynth_nufft: NUFFT complete in 0.00256 seconds.
INFO synthesis.get_rmsf_nufft: Fitting main lobe in each RMSF spectrum.
INFO rmsynth._run_rmsynth: RM-synthesis completed in 7.87ms.

[17]:

fdf_parameters

[17]:

shape: (1, 36)

fdf_error_mad	peak_pi_fit	peak_pi_error	peak_pi_fit_debias	peak_pi_fit_snr	peak_pi_fit_index	peak_rm_fit	peak_rm_fit_error	peak_q_fit	peak_u_fit	peak_pa_fit_deg	peak_pa_fit_deg_error	peak_pa0_fit_deg	peak_pa0_fit_deg_error	fit_function	lam_sq_0_m2	ref_freq_hz	fwhm_rmsf_radm2	phi_max_scale_radm2	fdf_error_noise	fdf_q_noise	fdf_u_noise	min_freq_hz	max_freq_hz	n_channels	median_d_freq_hz	frac_pol	frac_pol_error	sigma_add	sigma_add_minus	sigma_add_plus	mom0	mom0_debias	mom1_radm2	mom2_radm2	moment_threshold_snr
f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	str	f64	f64	f64	f64	f64	f64	f64	f64	f64	i64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64	f64
NaN	0.204894	0.01459	0.203696	14.043864	1265	89.364055	1.229083	-0.198363	-0.050445	97.134084	2.039887	34.397406	11.808332	"log"	0.082563	1.0433e9	34.522155	113.191071	0.01459	0.01461	0.014569	7.995e8	1.7995e9	54	8.2286e6	0.191109	0.013688	2.053029	1.584966	2.533219	0.623978	0.614144	115.404973	58.706994	5.0

[18]:

fdf_arrs

[18]:

shape: (2_001, 2)

phi_arr_radm2	fdf_dirty_complex_arr
f64	object
-336.725921	(-0.014737349935157283+0.009365828725447037j)
-336.389195	(-0.015268595888951995+0.009364196871523149j)
-336.052469	(-0.015792726350951483+0.009354890832883432j)
-335.715743	(-0.016309049885540203+0.009337566210316374j)
-335.379017	(-0.016816874886624596+0.009311858133831227j)
…	…
335.379017	(-0.02092803422270988+0.032044131457206144j)
335.715743	(-0.021480196465542785+0.031315291558238374j)
336.052469	(-0.022008208742584006+0.030587628749930814j)
336.389195	(-0.022511925225759278+0.029862374893173557j)
336.725921	(-0.022991231231907824+0.029140743757237148j)

[19]:

stokes_i_arrs

[19]:

shape: (54, 7)

freq_arr_hz	lambda_sq_arr_m2	stokes_i_model_arr	stokes_i_model_error	flag_arr	complex_pol_arr	complex_pol_error
f64	f64	f64	f64	bool	object	object
7.995e8	0.140606	1.275054	0.049528	true	(0.043383814093691325+0.0818397444768925j)	(0.07844614670963128+0.07849244435394578j)
8.0773e8	0.137756	1.266359	0.048256	true	(0.02395288434377864-0.16560746730739032j)	(0.07897182382436067+0.07921830430624699j)
8.1596e8	0.134992	1.257531	0.046509	true	(0.02458312861601038-0.027359994037060967j)	(0.07952608838960785+0.07952732881820707j)
8.2419e8	0.13231	1.249001	0.045496	true	(0.12765488594176416-0.05480416345919908j)	(0.08019889103739038+0.08008886093079004j)
8.3241e8	0.129707	1.240806	0.043504	true	(-0.0981080509377833-0.11250707602902542j)	(0.08066616385202692+0.08068926787314547j)
…	…	…	…	…	…	…
1.6555e9	0.032793	0.781783	0.048637	true	(0.5334475003168457+0.4421881290985238j)	(0.1321478635277816+0.13083749261472893j)
1.6915e9	0.031412	0.770445	0.049932	true	(0.4024660865524486+0.1253350788353451j)	(0.13238999908262558+0.13004902714864588j)
1.7275e9	0.030117	0.759815	0.0507	true	(0.5427981191609134-0.10499192143576891j)	(0.13650386780431548+0.13179737699333946j)
1.7635e9	0.0289	0.749301	0.051473	true	(0.4530268065449075-0.21893543202991322j)	(0.13703806558687562+0.1343024105153554j)
1.7995e9	0.027755	0.73941	0.052036	true	(0.5113292545828951-0.322674091057177j)	(0.13994840980790213+0.1371361599017985j)

[20]:

fig, ax = plt.subplots()
ax.plot(freq_hz, stokes_i_noiseless, label="Input model")
ax.plot(freq_hz, stokes_i_noisy, ".", label="Noisy data")
ax.plot(
    stokes_i_arrs["freq_arr_hz"],
    stokes_i_arrs["stokes_i_model_arr"],
    "k--",
    label="Fitted model",
)
ax.fill_between(
    stokes_i_arrs["freq_arr_hz"],
    stokes_i_arrs["stokes_i_model_arr"] - stokes_i_arrs["stokes_i_model_error"],
    stokes_i_arrs["stokes_i_model_arr"] + stokes_i_arrs["stokes_i_model_error"],
    alpha=0.3,
    color="k",
    label="Fitted model error",
)
ax.legend()
ax.set(
    xlabel=rf"$\nu$ / {u.Hz:latex_inline}",
    ylabel="Flux density",
    title="Stokes I",
)

[20]:

[Text(0.5, 0, '$\\nu$ / $\\mathrm{Hz}$'),
 Text(0, 0.5, 'Flux density'),
 Text(0.5, 1.0, 'Stokes I')]

[21]:

phi_arr_radm2 = fdf_arrs["phi_arr_radm2"].to_numpy()
fdf_dirty_arr = fdf_arrs["fdf_dirty_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()

ax.plot(
    phi_arr_radm2,
    fdf_dirty_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi_arr_radm2,
    fdf_dirty_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi_arr_radm2,
    np.abs(fdf_dirty_arr),
    color="k",
    label="Polarized intensity",
)

ax.errorbar(
    fdf_parameters["peak_rm_fit"],
    fdf_parameters["peak_pi_fit"],
    xerr=fdf_parameters["peak_rm_fit_error"],
    yerr=fdf_parameters["peak_pi_error"],
    fmt="o",
    lw=1,
    color="red",
    mfc="none",
    label="Fitted peak",
)

ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="Flux density",
    title="Dirty FDF",
)
ax.legend()

[21]:

<matplotlib.legend.Legend at 0x7d13a9166060>

[22]:

phi2_arr_radm2 = rmsf_arrs["phi2_arr_radm2"].to_numpy()
rmsf_arr = rmsf_arrs["rmsf_complex_arr"].to_numpy().astype(complex)

fig, ax = plt.subplots()
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.real,
    color="tab:red",
    label="Stokes Q",
)
ax.plot(
    phi2_arr_radm2,
    rmsf_arr.imag,
    color="tab:blue",
    label="Stokes U",
)
ax.plot(
    phi2_arr_radm2,
    np.abs(rmsf_arr),
    color="k",
    label="Polarized intensity",
)
ax.legend()
ax.set(
    xlabel=rf"$\phi$ / {u.rad / u.m**2:latex_inline}",
    ylabel="RMSF",
    title="RMSF",
)

[22]:

[Text(0.5, 0, '$\\phi$ / $\\mathrm{rad\\,m^{-2}}$'),
 Text(0, 0.5, 'RMSF'),
 Text(0.5, 1.0, 'RMSF')]