Electron-scale measurements of magnetic reconnection in space

1 MMS Science event: magnetopause reconnection at the electron diffusion region

Obtain and analyze plasma data, including spectra, for the first MMS Science event (Burch et al., Science 2016, link here) showing magnetopause reconnection at the electron diffusion region on 2016/10/16 13:07:02.2 UT. Start by using Hwk06_mpause_RX.pro, provided. The objective of this exercise is to introduce plasma distributions from MMS (electrons and ions) and create burst spectra. In the process also create plasma moments and the plasma current and plot these along E and B. You are requested to create Figure 2 of Burch et al. plus the 3 bottom (electron) spectrograms from Figure 3 of the same paper (panels 3G, 3H and 3I). Your figure should look like the one shown on the next page, but using burst mode data in order to make it look like in Burch et al. You are requested to plot this figure at 2 time scales: the overview (2 min) timescale as in Fig. 2 of Burch et al. [‘13:05:30’,’13:07:30’] and the zoom-in (3 sec) timescale as in Fig. 3 of Burch et al. [‘13:07:00.5’,‘13:07:03.5’]. Notice that the clean and fast m’pause crossing was at 13:05:40UT, and this is used to determine N. In Fig. 2K of Burch et al. (the right-hand side of Fig. 2, the pictorial view of the MMS trajectory for the 2min interval) this initial m’pause crossing was near the start of the trajectory. The trajectory crosses the X-point at the 3 seconds of the zoom-in interval.

Code 
dir = "docs/courses/epss261/homework"
if isdir(dir)
    cd(dir)
    Pkg.activate(".")
    Pkg.resolve()
    Pkg.instantiate()
end
Code 
using Dates
using Speasy
using DimensionalData
using SPEDAS
using SPEDAS.MMS
# using GLMakie
using Unitful
using CairoMakie, SpacePhysicsMakie

update_theme!(colormap=:rainbow1)
SpacePhysicsMakie.DEFAULTS.add_title = true;
Precompiling packages...
  10374.3 ms  ✓ SpacePhysicsMakie
  1 dependency successfully precompiled in 12 seconds. 297 already precompiled.

1.1 Overview plot

1.1.1 Minimum Variance Analysis

Code 
trange = ("2015-10-16T13:05:35", "2015-10-16T13:07:25")
tr_mpause = DateTime.(("2015-10-16T13:05:40", "2015-10-16T13:06:09"))

probe = 2

tvars = [
    "cda/MMS2_FGM_SRVY_L2/mms2_fgm_b_gse_srvy_l2_clean",
    "cda/MMS2_EDP_FAST_L2_DCE/mms2_edp_dce_gse_fast_l2",
]
B_gse, edp_dce_gse = get_data(tvars, trange) .|> DimArray
rotMat = mva_eigen(tclip(B_gse[:, 1:3], tr_mpause))
B_LMN = rotate(B_gse[:, 1:3], rotMat) |> set_coord("LMN") |> set_coord("Boundary-normal coordinates"; old_coords=["Geocentric Solar Ecliptic"])
E_LMN = rotate(edp_dce_gse, rotMat) |> set_coord("LMN")

rotMat
Can't get MMS2_FGM_SRVY_L2/mms2_fgm_b_gse_srvy_l2_clean without web service, switching to web service
LinearAlgebra.Eigen{Float32, Float32, StaticArraysCore.SMatrix{3, 3, Float32, 9}, StaticArraysCore.SVector{3, Float32}}
values:
3-element StaticArraysCore.SVector{3, Float32} with indices SOneTo(3):
 501.50854
  32.284042
  13.253316
vectors:
3×3 StaticArraysCore.SMatrix{3, 3, Float32, 9} with indices SOneTo(3)×SOneTo(3):
  0.368873   0.571603   0.732941
 -0.122903  -0.751631   0.648033
  0.921318  -0.329122  -0.207005

The direction obtained from the Minimum Variance Analysis (MVA) in our study closely aligns with the direction reported in the literature.

The (x, y, z) GSE components of the L, M, and N axes are L = (0.3665, –0.1201, 0.9226) GSE, M = (0.5694, –0.7553, –0.3245) GSE, and N = (0.7358, 0.6443, –0.2084) GSE

1.1.2 Plasma data

Code 
data_rate = "brst"
fpi_des_ds = FPIDataSet(; probe, data_rate, data_type="des")
fpi_dis_ds = FPIDataSet(; probe, data_rate, data_type="dis")

des_data = get_data(NamedTuple, fpi_des_ds, trange)
dis_data = get_data(NamedTuple, fpi_dis_ds, trange)

des_n = setmeta(DimArray(des_data.numberdensity)[:, 1], label="Ne")
dis_bulkv_lmn = rotate(DimArray(dis_data.bulkv_gse), rotMat) |> set_coord("LMN")
des_bulkv_lmn = rotate(DimArray(des_data.bulkv_gse), rotMat) |> set_coord("LMN")
dis_n = setmeta(DimArray(dis_data.numberdensity)[:, 1], label="Ni")
┌ 733-element DimArray{Float32, 1} mms2_dis_numberdensity_brst ┐
├──────────────────────────────────────────────────────────────┴───────── dims ┐
  ↓ Ti Sampled{UnixTimes.UnixTime} [2015-10-16T13:05:35.144103000, …, 2015-10-16T13:07:24.945661000] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 16 entries:
  "SCALETYP"        => "linear"
  "FILLVAL"         => Any[-1.0e31]
  "DEPEND_0"        => "Epoch"
  "FIELDNAM"        => "MMS2 FPI/DIS number density"
  "SI_CONVERSION"   => "1e6>m^-3"
  "VALIDMAX"        => Any[100000.0]
  :label            => "Ni"
  "DELTA_MINUS_VAR" => "mms2_dis_numberdensity_err_brst"
  "FORMAT"          => "E12.2"
  "VAR_TYPE"        => "data"
  "CATDESC"         => "MMS2 FPI/DIS ion number density during this burst"
  "LABLAXIS"        => "N"
  "DELTA_PLUS_VAR"  => "mms2_dis_numberdensity_err_brst"
  "DISPLAY_TYPE"    => "time_series"
  "VALIDMIN"        => Any[0.001]
  "UNITS"           => "cm^-3"
└──────────────────────────────────────────────────────────────────────────────┘
 ⋮      ⋱

1.1.3 Magnitudes of electron and ion convection velocities

Code 
Vi_perp_mag, Ve_perp_mag = let Vi = DimArray(dis_data.bulkv_gse), Ve = DimArray(des_data.bulkv_gse), B = B_gse
    tr = timerange(B)
    Vi_clip = tclip(Vi, tr)
    Ve_clip = tclip(Ve, tr)

    B_int_Vi = tinterp(B, Vi_clip)[:, 1:3]
    B_int_Ve = tinterp(B, Ve_clip)[:, 1:3]
    Vi_perp_mag = toproj(Vi_clip, B_int_Vi) |> tnorm
    Ve_perp_mag = toproj(Ve_clip, B_int_Ve) |> tnorm

    setmeta(Vi_perp_mag, ylabel="Viper_t"),
    setmeta(Ve_perp_mag, ylabel="Veper_t")
end
(Float32[16.370079, 39.08527, 18.779638, 36.385616, 28.834595, 31.232655, 21.860336, 25.042305, 19.904325, 3.4641075  …  54.85829, 52.707363, 43.587463, 37.501587, 32.176003, 29.350264, 22.573795, 23.164637, 26.775883, 33.856125], Float32[179.42912, 65.0223, 95.883736, 80.23359, 31.38685, 176.10417, 21.706322, 53.19559, 126.49119, 102.77766  …  42.916092, 72.561844, 39.365314, 40.899796, 34.341866, 35.257183, 31.910059, 45.77481, 69.92187, 93.317215])

1.1.4 Current density

Code 
J = begin
    Vi = dis_bulkv_lmn
    Ve = des_bulkv_lmn
    n = dis_n
    Ve_clip = tclip(Ve, timerange(Vi))
    Vi_interp = tinterp(Vi, Ve_clip)
    n_interp = tinterp(n, Ve_clip)
    J = mapslices(Vi_interp - Ve_clip; dims=Ti) do V_diff_i
        @. Unitful.q * V_diff_i * u"km/s" * n_interp * u"cm^-3" |> u"μA*m^-2"
    end
    setmeta(J, labels=["j_L", "j_M", "j_N"])
end
┌ 3661×3 DimArray{Unitful.Quantity{Float64, 𝐈 𝐋⁻², Unitful.FreeUnits{(μA, m⁻²), 𝐈 𝐋⁻², nothing}}, 2} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
  ↓ Ti Sampled{UnixTimes.UnixTime} [2015-10-16T13:05:35.144103000, …, 2015-10-16T13:07:24.945661000] ForwardOrdered Irregular Points,
  → AnonDim Sampled{Int64} [1, …, 3] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 19 entries:
  "UNITS"             => "km/s"
  "SCALETYP"          => "linear"
  "FILLVAL"           => Any[-1.0e31]
  "DEPEND_0"          => "Epoch"
  "FIELDNAM"          => "MMS2 FPI/DIS LMN bulk v"
  "SI_CONVERSION"     => "1.0e3>m s^-1"
  "VALIDMAX"          => Any[110000.0]
  "LABL_PTR_1"        => ["Vx_LMN", "Vy_LMN", "Vz_LMN"]
  "TENSOR_ORDER"      => Any[1]
  "COORDINATE_SYSTEM" => "LMN"
  "DELTA_MINUS_VAR"   => "mms2_dis_bulkv_err_brst"
  "FORMAT"            => "E12.2"
  "VAR_TYPE"          => "data"
  "CATDESC"           => "MMS2 FPI/DIS ion bulk-velocity LMN vector during this…
  :labels             => ["j_L", "j_M", "j_N"]
  "DELTA_PLUS_VAR"    => "mms2_dis_bulkv_err_brst"
  "DISPLAY_TYPE"      => "time_series"
  "VALIDMIN"          => Any[-110000.0]
  "REPRESENTATION_1"  => "mms2_dis_cartrep_brst"
└──────────────────────────────────────────────────────────────────────────────┘
 ⋮      ⋱

1.1.5 Para, perp and anti-parallel spectra of electrons

Code 
des_energyspectr_tvars = [
    "cda/MMS2_FPI_BRST_L2_DES-MOMS/mms$(probe)_des_energyspectr_par_$(data_rate)",
    "cda/MMS2_FPI_BRST_L2_DES-MOMS/mms$(probe)_des_energyspectr_perp_$(data_rate)",
    "cda/MMS2_FPI_BRST_L2_DES-MOMS/mms$(probe)_des_energyspectr_anti_$(data_rate)"
]

des_energyspectr = get_data(des_energyspectr_tvars, trange)
des_energyspectr = replace.(des_energyspectr, 0 => NaN)
des_energyspectr_ratio = des_energyspectr[1] ./ des_energyspectr[2]
des_energyspectr = setmeta.(des_energyspectr, colorrange=(1e5, 1e9))
des_energyspectr_ratio = setmeta(des_energyspectr_ratio,
    colorrange=(1e-1, 1e1), title="Ratio (Para/Perp)"
)
SpeasyVariable{Float64, 2}: mms2_des_energyspectr_par_brst
  Time Range: 2015-10-16T13:05:35.024103000 to 2015-10-16T13:07:24.975661000
  Units: 
  Size: (3666, 32)
  Memory Usage: 1.821 MiB
  Metadata:
    VAR_NOTES: Counts, summed within 30 degrees parallel bentPipe magnetic field.
    SCALETYP: log
    FILLVAL: Any[-9.999999848243207e30]
    DEPEND_0: Epoch
    FIELDNAM: MMS2 FPI/DES energySpectr_par
    VALIDMAX: Any[1.0000000150474662e30]
    colorrange: (0.1, 10.0)
    DEPEND_1: mms2_des_energy_brst
    FORMAT: E12.2
    VAR_TYPE: data
    CATDESC: MMS2 FPI/DES electron energy parallel spectrum 30 degrees parallel to B during this burst
    LABLAXIS: FPI1/DES EnSpectr, Parallel
    title: Ratio (Para/Perp)
    DISPLAY_TYPE: spectrogram
    VALIDMIN: Any[-1.0000000150474662e30]

1.1.6 Summary data plot

In 1 paragraph (10 lines) explain what each panel represents for this reconnection interval.

From top to bottom, the panels represent the following: 1. the magnetic field vectors in LMN coordinate system; 2. energy times spectrogram of ion energy flux; 3. energy times spectrogram of electron energy flux; 4. plasma density (ion and electron); 5. ion flow velocity vectors in LMN coordinate system; 6. magnitudes of electron and ion convection velocities; 7. current density; 8. electron parallel and perpendicular temperatures; 9. electric field vectors in LMN coordinate system; 10-12. electron spectrograms (parallel, perpendicular, anti-parallel); 13. electron spectrogram ratio (para/perp).

Code 
tvars2plot = [
    B_LMN,
    [dis_n, des_n],
    dis_data.energyspectr_omni,
    des_data.energyspectr_omni,
    dis_bulkv_lmn,
    [Vi_perp_mag, Ve_perp_mag],
    E_LMN,
    des_energyspectr...,
    des_energyspectr_ratio
]
faxes = tplot(tvars2plot)
tlines!(faxes, "2015-10-16T13:05:52")
tlines!(faxes, ["2015-10-16T13:05:44", "2015-10-16T13:06:00"])
tlines!(faxes, ["2015-10-16T13:06:46", "2015-10-16T13:07:12"])
ylims!.(faxes.axes[end-3:end], 1e1, 1e3)
faxes

1.2 Zoom-in electron spectrograms

Explain in 1 paragraph what the perpendicular electrons show.

Zoom in on the electron spectrograms

Code 
tBurchFig4 = ("2015-10-16T13:07:00.5", "2015-10-16T13:07:03.5")
tlims!(tBurchFig4)

The electron flow speed perpendicular to the magnetic field significantly exceeds the ion flow speed in the vicinity of the X-line, leading to a much stronger current near the X-line compared to the exhaust region. In contrast, the current closely follows the perpendicular ion speed within the magnetosheath and exhaust regions.

Furthermore, as illustrated in the figure below, reconnection dissipation is driven by the intense negative \(J_M\) current and negative \(E_M\) electric field, both perpendicular to the magnetic field \(B\) in the dissipation region. This condition preferentially accelerates electrons in the perpendicular direction, thereby reducing the ratio of electron parallel-to-perpendicular temperature. This effect is clearly demonstrated in the bottom panel of the figure above (compared to the magnetosheath and exhaust regions)

1.3 Dissipation quantity (J · E’)

Code 
B_interp = tinterp(B_LMN, J)
E_LMN_interp = tinterp(E_LMN, J) * u"mV/m"
J_dot_E = setmeta(tdot(J, E_LMN_interp), labels="j ⋅ E")
E_parp = setmeta(tsproj(E_LMN_interp, B_interp), ylabel="E∥ (mV/m)")
J_parp = tsproj(J, B_interp)
J_dot_E_para = setmeta(tdot(J_parp, E_parp), labels="j∥ ⋅ E∥")
J_dot_E_perp = setmeta(tdot(toproj(J, B_interp), toproj(E_LMN_interp, B_interp)), labels="j⊥ ⋅ E⊥")
┌ 3661-element DimArray{Unitful.Quantity{Float64, 𝐌 𝐋⁻¹ 𝐓⁻³, Unitful.FreeUnits{(μA, m⁻³, mV), 𝐌 𝐋⁻¹ 𝐓⁻³, nothing}}, 1} ┐
├──────────────────────────────────────────────────────────────────────── dims ┤
  ↓ Ti Sampled{UnixTimes.UnixTime} [2015-10-16T13:05:35.144103000, …, 2015-10-16T13:07:24.945661000] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 1 entry:
  :labels => "j⊥ ⋅ E⊥"
└──────────────────────────────────────────────────────────────────────────────┘
 2015-10-16T13:05:35.144103000  -0.0151287 μA mV m⁻³
 2015-10-16T13:05:35.174103000  -0.00870461 μA mV m⁻³
 2015-10-16T13:05:35.204103000   0.0100068 μA mV m⁻³
 ⋮                              
 2015-10-16T13:07:24.885661000   0.00910275 μA mV m⁻³
 2015-10-16T13:07:24.915661000   0.0203746 μA mV m⁻³
 2015-10-16T13:07:24.945661000   0.0440117 μA mV m⁻³

J · E is the total energy transfer (energy conversion) rate, which is a key quantity in reconnection studies. Since reconnection is known to be a dissipative process that converts magnetic energy to heat and particle kinetic energy, the dissipation quantity could help identify a physically relevant, small-scale region surrounding the reconnection.

Code 
tvars2plot = [
    B_interp,
    J,
    E_LMN,
    E_parp,
    [J_dot_E, J_dot_E_para, J_dot_E_perp],
]
faxes = tplot(tvars2plot)
tlines!(faxes, "2015-10-16T13:05:52")
tlines!(faxes, ["2015-10-16T13:05:44", "2015-10-16T13:06:00"])
tlines!(faxes, ["2015-10-16T13:06:46", "2015-10-16T13:07:12"])
faxes

Zoom-in on the dissipation region

Code 
tlims!(tBurchFig4)