© The Authors 2024

1. Introduction

The degree of polarisation in spectral lines induced by magnetic fields in the solar atmosphere can be higher than 10% in regions with strong field (e.g. Martinez Pillet et al. 1990; Castellanos Durán et al. 2020), but more commonly it is much lower: linear polarisation induced by the Zeeman and Hanle effect in lines formed in the chromosphere tends to be smaller than 0.1%, especially in the quiet Sun where field strengths are low (e.g. Jurčák et al. 2018; Centeno et al. 2021). Observations of the full Stokes vector therefore require long integration times and/or spatial averaging over areas larger than the diffraction limit of the telescope in order to reach the signal-to-noise ratio (S/N) required to detect signals with a strength of < 0.1% of the local continuum intensity (e.g. Centeno et al. 2010).

A major goal of modern solar telescopes is to perform such imaging spectropolarimetry at the diffraction limit. This goal is largely reliant on the calibration accuracy of the instruments. For example, the CRISP instrument at the Swedish 1-m Solar Telescope (Scharmer et al. 2003) can currently reach a polarimetric accuracy of ∼10⁻³ and the DKIST telescope can be calibrated to an accuracy of 5 × 10⁻⁴ (Rimmele et al. 2020).

Because of the finite brightness of the Sun, the S/N that is required to match the polarimetric accuracy is not reached instantaneously. For example, it takes 2 s to reach a S/N of 5 × 10⁻⁴ at the diffraction limit in the continuum adjacent to the Ca II 854.2 nm line at a spectral resolution of λ/Δλ = 80 000 and a total transmission of the telescope-instrument-camera chain of 5%. Tuneable-filter-based instruments, such as CRISP or VTF at DKIST, require data acquisition typically at 10 to 20 wavelengths, depending on the width of the line and the desired wavelength coverage. This leads to a total acquisition time of roughly 20 to 40 s for a diffraction-limited full-Stokes line scan.

This time is large compared to the timescale on which the density, temperature, and velocity in the atmosphere change: the sound-crossing time for a diffraction-limited pixel of a 1 m telescope at 854 nm is 9 s assuming a sound speed of 7 km s⁻¹; for a 4 m telescope, this value is 2 s. The sound speed should be considered as a lower limit; in the chromosphere, the Alfvén speed is a better indicator of the evolution speed. This latter can be larger than 100 km s⁻¹, implying changes of the thermodynamic parameters in less than 0.1 s.

Therefore, using diffraction-limited instruments based on tunable filters to reach the required high S/N inevitably leads to observed Stokes profiles that average different atmospheric states in density, temperature, and velocity. Using integral field spectropolarimeters, such as microlensed hyperspectral imagers (van Noort et al. 2022) or image slicers (Dominguez-Tagle et al. 2022), mitigates – but does not fully eliminate – this effect, but at the cost of a strongly reduced field of view (FOV).

This mismatch of required integration times and evolution timescales poses a severe problem for inferring chromospheric magnetic fields, especially over the large FOVs of Fabry-Pérot-type instruments. The equations that describe radiative transfer are non-linear. An inversion of a time-averaged line profile does not, in general, yield an inferred atmosphere that represents the time-average of the underlying atmosphere (Carlsson & Stein 1995). Time-averaging of the Stokes vector might even result in an observed average Stokes vector that cannot be reproduced by any single atmosphere model.

The magnetic field in the chromosphere, where the magnetic pressure is higher than the gas pressure, tends to be relatively smooth and to vary slowly over space. Morosin et al. (2020) realised that this property can be exploited to improve the inference of magnetic fields from noisy data using the weak-field approximation (WFA). These authors extended the merit function used in the WFA-fit with a term that penalises differences in magnetic field in neighbouring pixels, and showed that this leads to substantially improved inference of the longitudinal magnetic field based on the Mg I 517.3 nm, Na I 589.6 nm, and Ca II 854.2 nm lines.

The magnetic field in the chromosphere generally only evolves over granular turnover timescales of several minutes (Kleint 2017; Ferrente et al. 2023), but the mass density, temperature, and velocity evolve with a timescale of seconds on spatial scales of a few tens of kilometers. We propose to exploit this temporal smoothness to improve magnetic field inference from observations taken at or close to the diffraction limit.

The general idea is to perform inversions where, for the same spatial pixel, observations at consecutive time steps are allowed to have arbitrary variations in the velocity, temperature, and other variables that are well-constrained by Stokes I at relatively low S/N. The magnetic field, which is constrained by Stokes Q, U, and V and requires a high S/N to be constrained, should vary as little as possible while staying consistent with the data, given the noise level.

This concept can be viewed as applying Tikhonov regularisation (Tikhonov & Arsenin 1977; Willoughby 1979) to the magnetic field in the time domain. In the case of the WFA, this can be implemented as a straightforward extension of the method presented by Morosin et al. (2020). In the case of inference methods based on explicit radiative transfer, whether in the Milne-Eddington approximation (ME), LTE, or full non-LTE, inference methods can be implemented as an extension of the method described in de la Cruz Rodríguez (2019). The WFA is linear, and requires only a single matrix inversion to yield a solution; while the rest of the aforementioned inference methods are non-linear and can be solved using the Levenberg-Marquardt algorithm (LM, Levenberg 1944; Marquardt 1963).

In the present paper, we provide implementations for both cases. As a proof of concept, we apply the temporally regularised WFA to observations of quiet Sun in the Ca II 854.2 nm line, and apply temporally regularised ME inversions to observations of quiet Sun in the Fe I 617.3 nm line. We chose to illustrate the non-linear case with the ME approximation applied to a photospheric line because it is easy to implement and can be quickly inverted; the formalism also remains identical for LTE and non-LTE inversions of chromospheric lines.

The observations that we use here are obtained with a Fabry-Perot tunable-filter instrument. The individual line positions contained in the data are not taken simultaneously, but sequentially. We assume that one line scan comprises a time step that will be represented by one atmosphere model. As we mention above, this is an approximation, as the solar evolution timescale is smaller than the scan time, and this can translate into errors in the inferred model parameters (e.g. Felipe et al. 2018). We stress that temporal regularisation can also be used for integral-field spectropolarimeters. These record all wavelengths simultaneously, and the integration time used in a time step is mainly set by the desired S/N, and this time is significantly shorter than would be needed by a tunable-filter instrument.

2. Observations

We use two datasets acquired with the CRISP instrument (see Scharmer 2006; Scharmer et al. 2019; de Wijn et al. 2021) at the Swedish 1-m Solar Telescope. The datasets were reduced with the SSTRED pipeline (de la Cruz Rodríguez et al. 2015; Löfdahl et al. 2021), and processed with the MOMFBD algorithm to minimise atmospheric distortions (Löfdahl 2002; Van Noort et al. 2005). The polarimetric calibration was individually performed for each pixel as described in van Noort & Rouppe van der Voort (2008).

The first dataset, recorded on 14 July 2020 starting at 08:40 UT at disc centre, consists of a four-line program with data acquired in the Ca II 854.2 nm, Fe I 617.3 nm, Hα and Mg I 517.3 nm lines. We only used the 854.2 nm data in this study. The line was sampled between ±30 pm from line centre in regular steps of 7.5 pm, with a total integration time per wavelength of 0.91 s. The overall cadence – including all four lines – is 34.2 s, whereas the total acquisition time for one 854.2 nm scan was approximately 13 s. The upper row in Fig. 1 illustrates the four Stokes parameters in the first time step at +15 pm from line centre. The target is a quiet-Sun region that harbours two small network patches with opposite polarity. We estimated the noise level based on the outermost wavelength points. The standard deviations of the noise were found to be σ_Q, U, V = (3.9,  2.6,  2.7)×10⁻³.

Fig. 1. Observations used in this study. Top: Stokes parameters in the Ca II 854.2 nm line at +15 pm from line centre. Bottom: Stokes parameters in the Fe I 617.3 nm line at −17.5 (I) and −7.0 pm (Q, U, V) from line centre. Stokes Q, U, and V have been normalised by the spatially averaged Stokes I intensity at the same wavelength.

The second dataset, recorded on 26 May 2021 starting at 09:48 UT on coordinates (X, Y) = (698, 387) (μ = 0.54), consists of a patch of quiet Sun close to active region AR 12824. We recorded data in the Ca II 854.2 nm and Fe I 617.3 nm lines, but we only used the 617.3 nm data in this study. The Fe I 617.3 nm line was sampled at Δλ = [ − 17.5,   − 10.5,   − 7,   − 3.5,  0,  3.5,  7,  10.5,  14,  17.5] pm from line centre, with a total integration time per wavelength of 0.28 s. The selected FOV is illustrated in the bottom row of Fig. 1 in all four Stokes parameters at −7 pm from line centre. The background noise as estimated from the continuum point has standard deviations of σ_Q, U, V = (1.8,  2.3,  2.5)×10⁻³.

3. Method

3.1. The linear case: Application to the weak-field approximation

When the dependence of the observable on the model parameters is linear, spatial and temporal regularisation can be trivially imposed. The weak-field approximation (Landi Degl’Innocenti & Landi Degl’Innocenti 1973; Jefferies et al. 1989) allows us to write analytical expressions for the Q, U, and V profiles that are linear in $B_{\perp }^2$ and B_∥.

We extend the implementation of Morosin et al. (2020) by also adding regularisation in the temporal dimension. In order to derive general expressions for the linear case, let us assume that we can express our synthetic prediction s_i of the ith data point in a pixel as s_i = c_iP, where c_i is a quantity that does not depend on the model parameters (but that can change for each data point) and P is a model parameter. For a given pixel at coordinates t, y, x, we can write the regularised merit function χ² as:

$$\begin{array}{cc}\hfill {\chi }^2=& \frac{1}{N_{\mathrm{data}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{data}}}}{\left(\frac{o_{t,y,x}^i-c_i{P}_{t,y,x}}{{\sigma }_i}\right)}^2+\alpha [{(P_{t,y,x}-P_{t,y,x-1})}^2\hfill \\ \hfill & +{(P_{t,y,x}-P_{t,y,x+1})}^2+{(P_{t,y,x}-P_{t,y-1,x})}^2+{(P_{t,y,x}-P_{t,y+1,x})}^2]\hfill \\ \hfill & +\beta {(P_{t,y,x}-P_0)}^2+\gamma [{(P_{t,y,x}-P_{t-1,y,x})}^2+{(P_{t,y,x}-P_{t+1,y,x})}^2],\hfill \end{array}$$(1)

where o_i is the ith observed data point, α is a spatial regularisation weight, β is a low-norm regularisation weight that penalises deviations from a constant value P₀ and γ is a temporal regularisation weight. Taking the derivative of Eq. (1) with respect to all P_t, y, x, and equating it to zero to find the minimum, while absorbing the scaling factors N_data and σ_i into o_i and c_i, we can derive a sparse linear system of equations (AP = B), where the different temporal and spatial locations are coupled by the non-diagonal regularisation terms:

$$\begin{array}{cc}\hfill [{\displaystyle \sum _i}{c}_i^2+4\alpha +\beta ]P_{t,y,x}& -\alpha (P_{t,y-1,x}+P_{t,y+1,x}+P_{t,y,x-1}+P_{t,y,x+1})\hfill \\ \hfill & -\gamma (P_{t-1,y,x}+P_{t+1,y,x})={\displaystyle \sum _i}{c}_i{o}_{t,y,x}^i+\beta P_0.\hfill \end{array}$$(2)

Equation (2) couples the parameters of the different pixels (P_t, y, x) through the regularisation weights (α, β, γ) on the left-hand side. Obtaining a solution requires that we solve a global problem for all pixels simultaneously. If the regularisation weights are zero, the A matrix in Eq. (2) is diagonal and we recover the unconstrained algorithm. Figure 2 shows the structure of A for a small FOV of 4 × 4 pixels and three time steps. Compared to Fig. 1 of Morosin et al. (2020), this matrix has two extra off-diagonal bands originating from the temporal regularisation terms.

Fig. 2. Regularisation matrix for the linear case. Each row corresponds to the regularisation function for one pixel in space and time. For displaying purposes, we assume an observation with dimensions nt = 3, ny = 4, nx = 4 and regularisation weights α = 1, β = 1, and γ = 2.5. The outermost bands (darker blue) correspond to the temporal regularisation terms, whereas the inner light-blue terms originate from the spatial regularisation.

If we replace P_t, y, x with $B_{\Vert }^{t,y,x}$ and ∑_ic_i with ∑_iC∂I_i/∂λ_i, we directly obtain the expression for the fully regularised WFA, with $C=4.67\times {10}^{-12}{\lambda }_0^2{g}_{\text{eff}}$, λ₀ the central wavelength of the line in nm, and g_eff the effective Landé factor of the transition.

Figure 3 illustrates a comparison of four reconstructions of B_∥ from the Ca II 854.2 nm dataset using the unconstrained WFA algorithm, the spatially regularised method of Morosin et al. (2020), and two variants of our method that include temporal regularisation alone as well as both temporal and spatial regularisation. The reconstruction of the unconstrained algorithm has a very noisy background, and besides the magnetic field in the strong network patches, small-scale loop-like features can be discerned in the FOV. The inclusion of spatial regularisation greatly decreases the noise in the background, making the small-scale loop-like features much more visible compared to the unconstrained case. Temporal regularisation alone also decreases the noise, perhaps yielding a slightly sharper model than the spatially regularised one, while also decreasing the temporal fluctuations of the background noise along the time series. The combined action of temporal and spatial regularisation further decreases the noise compared to the previous cases and yields a model with the highest S/N. The standard deviation of the noise in the magnetic field, as measured in the lower-right corner of the image, is reduced from 18 G in the unconstrained case to 9 G in the fully regularised case. The improvement induced by spatial regularisation is particularly obvious in the movie showing the entire time series. The upper panel in Fig. 4 further illustrates the reduction of the noise in the regularised reconstructions along a slice through the FOV.

Fig. 3. Reconstructions of the line-of-sight component of the magnetic field (B∥), inferred by applying the weak-field approximation to the Ca II 854.2 nm data. Upper-left panel: reconstruction performed with the unconstrained algorithm. Lower-left panel: reconstruction performed with a temporally regularised algorithm. Lower-right panel: reconstruction performed with a spatially regularised algorithm. Upper-right panel: reconstruction performed with full regularisation in space and time. A movie with the time series is available https://www.aanda.org/10.1051/0004-6361/202348810/olm. The red tick marks indicate the location of the cut shown in Fig. 4.

Fig. 4. Vertical cuts of the reconstructed magnetic field. The upper panel shows a reconstruction of B∥ with the weak-field approximation from the 854.2 nm dataset along a cut indicated with red markers in Fig. 3. The black curve shows the unconstrained reconstruction, and the red curve shows the reconstruction with both spatial and temporal regularisation. A similar plot is shown in the lower panel for the reconstruction of B⊥ from the ME inversion of the 617.3 nm dataset presented in Fig. 7.

3.2. The non-linear case: Application to Milne-Eddington inversions

The LM algorithm is one of the most efficient methods for reconstructing the parameters of non-linear models from observations. Regularised LM implementations are used in different fields of astrophysics in order to constrain the parameters of the model under consideration (e.g. Piskunov & Kochukhov 2002; de la Cruz Rodríguez et al. 2019). The regularisation is represented by a set of N_pen penalty functions r_n(p) that are squared in order to have ℓ − 2 norm regularisation. Following the notation of de la Cruz Rodríguez et al. (2019), the merit function χ² can generally be expressed as a function of the model parameter vector p and the data points x:

$$\begin{array}{c}\hfill {\chi }^2(\mathit{p},\mathit{x})=\frac{1}{N_{\mathrm{data}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{data}}}}{\left[\frac{o_k-s(\mathit{p},x_k)}{{\sigma }_k}\right]}^2+{\displaystyle \sum _{n=1}^{N_{\mathrm{pen}}}}{\alpha }_n{r}_n{(\mathit{p})}^2,\end{array}$$(3)

where the weights α_n regulate the influence of the regularisation terms in the merit function. We note that, unlike in the linear case, the penalty functions are not independently defined for each pixel (see below).

The model corrections predicted by the regularised LM algorithm can be derived by linearising Eq. (3) and taking the derivative with respect to the model parameters. The correction to the model parameters (Δp), in vector form, is given by a linear system of equations:

$$\begin{array}{c}\hfill (\mathbf{J}·{\mathbf{J}}^T+\mathbf{L}·{\mathbf{L}}^T)\mathrm{\Delta }\mathit{p}=\mathbf{J}·(\mathit{o}-\mathit{s})-\mathbf{L}·\mathit{r},\end{array}$$(4)

where J is the Jacobian matrix of the model and L is the Jacobian matrix of the penalty functions. In this expression, all constants and normalising factors are implicitly contained in the corresponding vectors. For any given parameter p_t, y, x, the regularisation functions are defined as:

$$\begin{array}{cc}\hfill {\displaystyle \sum {\mathit{r}}_{t,y,x}^2}=& {\alpha }_t{(p_{t,y,x}-p_{t-1,y,x})}^2+{\alpha }_s{(p_{t,y,x}-p_{t,y-1,x})}^2\hfill \\ \hfill & +{\alpha }_s{(p_{t,y,x}-p_{t,y,x-1})}^2+{\alpha }_l{(p_{t,y,x}-p_0)}^2,\hfill \end{array}$$(5)

where we have included spatial, temporal, and low-norm regularisation terms.

If the penalty functions have a linear dependence on the model parameters, as in the ones used in this study, we can express them as r = Lp + c, where the coupling is contained in the Jacobian matrix L. Figure 5 illustrates the structure of L for a problem with n_t = 3, n_y = 4, n_x = 4, and n_par = 2. The maximum number of penalty functions should be n_pen < 4n_xn_yn_t, as some functions are not defined at the edges of the problem. Figure 6 illustrates the approximate Hessian matrix of the regularisation functions L ⋅ L^T. The size of the approximate Hessian matrix is set by the total number of free parameters of the problem, which is also the size of one row of the L^T matrix (4 × 4 × 3 × 2 = 96). In the latter, the outer dark blue bands originate from the temporal regularisation, whereas the inner light-blue bands contain the spatial coupling terms. For n_par = 1, this matrix is identical in form to that shown in Fig. 2 for the linear case.

Fig. 5. Transpose of the Jacobian matrix of the penalty functions for the non-linear case with nt = 3, ny = 4, nx = 4, npar = 2, and regularisation weights α = 1, β = 1, and γ = 2.5. Each row corresponds to one penalty function and each column indicates the free parameters that are affected.

Fig. 6. Hessian approximation of the regularisation functions (LTL). This matrix was computed using L from Fig. 5.

In this case, we are only including penalty functions that compare a parameter value with its neighbouring values at t − 1, y − 1, and x − 1. However, because the product L ⋅ L^T is similar to a correlation of each penalty function with all others, the resulting coupling matrix in the left-hand side has identical structure to that in the linear case where we also include explicit comparisons with the values at t + 1, y + 1, and x + 1.

In a previous study, de la Cruz Rodríguez (2019) implemented a multi-resolution inversion code that included spatial regularisation and allowed us to deal with the effects of the telescope point spread function in ME inversions (PyMilne¹). We extended this code to include temporal regularisation terms. With these changes, the code can now invert a time series of maps as a global problem. We inverted a time-series of 14 line scans acquired in the Fe I 617.3 nm line in a quiet-Sun patch. These observations are not ideal because they also include scans in the Ca II 854.2 nm line, inducing a time gap between consecutive Fe I scans. Nevertheless, the cadence is sufficiently high to illustrate the advantages of the method. In the following, we focus on the reconstruction of B_⊥, which is a more challenging problem than reconstructing B_∥.

We performed reconstructions with the unconstrained algorithm, including spatial regularisation only, temporal regularisation only, and a case including low-norm, spatial, and temporal regularisation altogether. Regularisation was only imposed in the three components of the magnetic field, leaving all other parameters unconstrained. Therefore, the thermodynamical quantities of the model are allowed to change freely in time and space.

Figure 7 compares reconstructions of the time series of B_⊥ from the unconstrained algorithm and the three regularised cases. The improvements of the fully regularised reconstruction are threefold:

Fig. 7. Reconstruction of the transverse component of the magnetic field vector from ME inversions of the Fe I 617.3 nm. Left: unconstrained reconstruction where each pixel has been inverted six times with five different random initialisations. Middle-left: reconstruction with spatial regularisation, where all pixels have only been inverted once. Middle-right: reconstruction with temporal regularisation. Right: reconstruction including temporal, spatial, and low-norm regularisation terms. A movie with the time-evolution is available https://www.aanda.org/10.1051/0004-6361/202348810/olm.

• So-called inversion noise – caused by pixels where the algorithm did not converge to the global minimum – is absent. This effect is clearest in the patches of strong magnetic field, where randomly located black pixels are present in the unconstrained reconstruction. Examples of inversion noise are also visible in the lower panel of Fig. 4 at y ≈ 3.7″ and y ≈ 4.35″. Inversion noise can be removed from the unconstrained method by initialising the inversion several times with different starting model parameters or by performing a two-cycle inversion with a smoothing step in between (as suggested by Milić & van Noort 2018). However, the noise is naturally addressed in all the regularised cases, where the coupling between pixels in space and/or time removes it entirely.

• There is a reduction of the background noise bias level. The presence of noise in the observations can manifest as a non-zero background bias-level for parameters that are defined positive (see Martínez González et al. 2012). We find that all forms of regularisation decrease this bias, but the low-norm regularisation term (with P₀ = 0) and the temporal regularisation are particularly effective in this regard. The spatially regularised time series shows a stronger frame-to-frame fluctuation in the background level, which seems to be induced by varying atmospheric seeing conditions along the series. This effect is also visible in the lower panel of Fig. 4, where the red curve is regularly predicting much lower values of B_⊥ in regions located between stronger network patches.

• The time coherence of weak magnetic structures is enhanced. Structures with signal barely above the noise appear more persistent in time with temporal regularisation. Without this regularisation, they appear intermittently, only in those time steps where the noise happens to be low enough for the inversion to pick up the structure.

However, despite these improvements, the evolution of magnetic fields in the photosphere is anchored to the forces originating from gas-pressure gradients (plasma β ≫ 1), and our assumption that the magnetic field is slowly evolving relative to the other parameters is not optimal. Compared to the chromospheric case presented in Sect. 3.1, we had to limit the amount of temporal regularisation to perform these inversions.

4. Discussion and conclusions

We present a method to include Tikhonov regularisation in the time dimension for linear and non-linear fits of model parameters. We show that this method is suitable for improving determinations of the magnetic field in the solar chromosphere.

Compared to other techniques that operate on the data to decrease noise, such as temporal and spatial rebinning, principal component analysis filtering (see e.g. Martínez González et al. 2008), Fourier filtering, or neural-network-based noise reduction (Díaz Baso et al. 2019), regularisation leaves the data untouched and operates only on the model parameters. This is an important advantage because filtering and/or averaging the data likely affects the reconstruction of all parameters, whether that is desired or not.

Temporal regularisation addresses the different timescales on which the magnetic field and the thermodynamical parameters evolve in the chromosphere as well as the different S/N required to constrain them. Because the thermodynamic parameters evolve relatively quickly but require a low S/N, one can use zero (as we do here) or weak regularisation. The magnetic field evolves more slowly but requires high S/N, and is therefore regularised more strongly.

Choosing adequate regularisation weights is crucial to obtain results that are not overly smooth. The weights must be chosen in such a way that the final quality of the fit is not affected; a threshold value should be chosen above which the quality of the fits starts to be negatively affected by the limitations imposed by the regularisation terms (Kochukhov 2017). For observations with low S/N, this might mean that the regularisation is so strong that it makes the magnetic field almost constant in time. We emphasise that with properly chosen weights, regularisation does not guarantee that the magnetic field will be constant over time; if the S/N is sufficiently high, it still allows sharp gradients in time. However, in cases of intermediate S/N, regularisation will produce smoother solutions if strong gradients are not required to minimise χ². Regularisation has an advantage over smoothing the model parameters after an unconstrained inversion, because smoothing washes out variations regardless of the impact that they have on χ², but regularisation maintains sharp gradients in locations where they have a significant impact in χ².

Most observational programs with Fabry-Perot instruments that observe multiple lines are carried out by scanning once through each line and repeating this cycle. Long integration times are needed for chromospheric lines in order to reach sufficient S/N to measure polarisation. As mentioned in the introduction, this takes about 20 s, and the observed line profile is not produced by a single atmosphere as the Sun has evolved over this time.

Temporal regularisation allows observing strategies that mitigate this problem. Instead of a single slow scan, one can perform multiple fast scans that each have sufficient S/N to constrain the thermodynamic parameters. These profiles now suffer less from solar evolution and therefore lead to improved reconstructions of the thermodynamic parameters. Temporal regularisation of the magnetic field couples the scans together and increases the “effective” S/N sufficiently to determine the field. These fast scans can be interleaved with observations in other lines (which themselves can consist of multiple scans) as needed, as temporal regularisation does not require a constant cadence. If needed, the temporal regularisation weight can be made smaller for scans that are not performed back-to-back.

Another obvious application is to data taken with integral-field spectropolarimeters. Here the assumptions behind temporal regularisation hold much better, because all spatial and spectral datapoints are taken simultaneously, and the cadence is typically high (Rouppe van der Voort et al. 2023).

Data taken at the diffraction limit of a telescope have a limited S/N if the exposure time is kept short enough to prevent solar evolution from blurring the images. Spatio-temporal regularisation of the magnetic field is a way of increasing the effective S/N of the data while keeping exposure times short. This technique will be invaluable when considering data taken with large-aperture telescopes, such as DKIST (Rimmele et al. 2020) and the planned EST (Quintero Noda et al. 2022).

The version of pyMilne that includes temporal regularisation is available in the official repository¹. The regularised WFA code used here is also publicly available².

Movies

Movie 1 associated with Fig. 3 (WFA) Access here

Movie 2 associated with Fig. 7 (ME_inversion) Access here

Acknowledgments

The authors of this paper are thankful to the referee for his careful evaluation of the manuscript. The Institute for Solar Physics is supported by a grant for research infrastructures of national importance from the Swedish Research Council (registration number 2021-00169). The Swedish 1-m Solar Telescope is operated on the island of La Palma by the Institute for Solar Physics of Stockholm University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. J.L. acknowledges financial support from the Swedish Research council (VR, project number 2022-03535). This project has been funded by the European Union through the European Research Council (ERC) under the Horizon Europe program (MAGHEAT, grant agreement 101088184).

References

All Figures

	Fig. 1. Observations used in this study. Top: Stokes parameters in the Ca II 854.2 nm line at +15 pm from line centre. Bottom: Stokes parameters in the Fe I 617.3 nm line at −17.5 (I) and −7.0 pm (Q, U, V) from line centre. Stokes Q, U, and V have been normalised by the spatially averaged Stokes I intensity at the same wavelength.
In the text

Fig. 2. Regularisation matrix for the linear case. Each row corresponds to the regularisation function for one pixel in space and time. For displaying purposes, we assume an observation with dimensions nt = 3, ny = 4, nx = 4 and regularisation weights α = 1, β = 1, and γ = 2.5. The outermost bands (darker blue) correspond to the temporal regularisation terms, whereas the inner light-blue terms originate from the spatial regularisation.

In the text

Fig. 3. Reconstructions of the line-of-sight component of the magnetic field (B∥), inferred by applying the weak-field approximation to the Ca II 854.2 nm data. Upper-left panel: reconstruction performed with the unconstrained algorithm. Lower-left panel: reconstruction performed with a temporally regularised algorithm. Lower-right panel: reconstruction performed with a spatially regularised algorithm. Upper-right panel: reconstruction performed with full regularisation in space and time. A movie with the time series is available https://www.aanda.org/10.1051/0004-6361/202348810/olm. The red tick marks indicate the location of the cut shown in Fig. 4.

In the text

Fig. 4. Vertical cuts of the reconstructed magnetic field. The upper panel shows a reconstruction of B∥ with the weak-field approximation from the 854.2 nm dataset along a cut indicated with red markers in Fig. 3. The black curve shows the unconstrained reconstruction, and the red curve shows the reconstruction with both spatial and temporal regularisation. A similar plot is shown in the lower panel for the reconstruction of B⊥ from the ME inversion of the 617.3 nm dataset presented in Fig. 7.

In the text

	Fig. 5. Transpose of the Jacobian matrix of the penalty functions for the non-linear case with nt = 3, ny = 4, nx = 4, npar = 2, and regularisation weights α = 1, β = 1, and γ = 2.5. Each row corresponds to one penalty function and each column indicates the free parameters that are affected.
In the text

	Fig. 6. Hessian approximation of the regularisation functions (LTL). This matrix was computed using L from Fig. 5.
In the text

Fig. 7. Reconstruction of the transverse component of the magnetic field vector from ME inversions of the Fe I 617.3 nm. Left: unconstrained reconstruction where each pixel has been inverted six times with five different random initialisations. Middle-left: reconstruction with spatial regularisation, where all pixels have only been inverted once. Middle-right: reconstruction with temporal regularisation. Right: reconstruction including temporal, spatial, and low-norm regularisation terms. A movie with the time-evolution is available https://www.aanda.org/10.1051/0004-6361/202348810/olm.

In the text