© The Authors 2026

1. Introduction

Through the widely employed techniques of spectroscopy and spectropolarimetry (which rely on the Doppler effect), it is only possible to infer the line-of-sight velocity of solar plasma. For simplicity, we refer to this velocity component, using Cartesian coordinates, as v_z at the disk center. However, many interesting phenomena that occur in the solar atmosphere involve the presence of large velocities in the plane perpendicular to the observer’s line of sight. We refer to these as v_x and v_y components of the velocity field. Some of these phenomena include flows at a sub-granular scale (Steiner et al. 2010), the presence of swirling motions (Bonet et al. 2010), diverging flows from magnetic reconnection events (Thaler & Borrero 2023), and magnetic tornadoes (Wedemeyer-Böhm et al. 2012). Apart from these phenomena, determining horizontal velocities (v_x, v_y) is extremely important because they can be employed to estimate the flux of magnetic energy and helicity from the photosphere into the corona (Welsch et al. 2007; Fu & Welsch 2014; Kazachenko et al. 2015; Tilipman et al. 2023) as well as to improve the inference of electric currents in the solar atmosphere (Pastor Yabar et al. 2021; Borrero & Pastor Yabar 2023).

The most widely used technique to determine velocities in the solar atmosphere in the plane perpendicular to the line of sight is the local correlation tracking technique (LCT; November & Simon 1988; Simon et al. 1988). This is an optical technique that uses time-resolved monochromatic intensity images, typically at a continuum wavelength where no spectral lines are present. The LCT has been successfully employed to determine motions on the solar surface in a variety of photospheric structures, such as granulation (Simon et al. 1995), pores (Hirzberger et al. 2002), sunspot penumbra (Denker 1998; Sobotka & Sütterlin 2001), sunspot moats (Vargas Domínguez et al. 2008), and vortices (Bonet et al. 2010), among others. Newer techniques include the use of deep learning neural networks (Asensio Ramos et al. 2017) and a local correlation tracking technique based on Fourier analysis (FLCT), which has been applied to time-resolved images of the continuum intensity (Fisher & Welsch 2008) as well as to time-resolved maps of the line-of-sight component of the magnetic field B_z (Welsch et al. 2004; DeForest et al. 2007; Welsch 2015; Kostić et al. 2026).

Assuming that the horizontal velocity inferred using LCT, referred to as u_h, is equivalent to real plasma velocities u_h = v_h (Kusano et al. 2002) or related to them via u_h = v_h − (v_z/B_z)B_h, for example, as proposed by Démoulin & Berger (2003), the vertical component (i.e., z-component) of the induction equation in ideal magnetohydrodynamics (MHD) can be incorporated to determine the full velocity vector v using different techniques (Welsch et al. 2004; Longcope 2004; Schuck 2008). These are referred to as inductive techniques, and they assume that the time derivative of B_z is known (i.e., obtained from a time series of line-of-sight magnetograms).

All the aforementioned methods share a common feature: They all retrieve the full velocity vector¹v in a (x, y) plane at a fixed height (or single z height) on the solar atmosphere: v(x, y). In this paper, we present an extension and modification of these methods where all three components of the induction equation are employed, enabling the inference of the horizontal components of the velocity at multiple heights in the solar atmosphere: v_h(x, y, z).

This new method presumes that the full magnetic field vector B(x, y, z, t) and the line-of-sight component of the plasma velocity v_z(x, y, z, t) are known. Nowadays, these physical parameters can be inferred, both in the photosphere (Borrero et al. 2025) and chromosphere (Kaithakkal et al. 2023), through advanced inversion techniques for the polarized radiative transfer equation (Pastor Yabar et al. 2019; Borrero et al. 2021) when applied to observations of the Stokes vector in photospheric and chromospheric spectral lines as a function of time, I(x, y, λ, t), and carried out with modern spectropolarimeters (Dominguez-Tagle et al. 2022; van Noort et al. 2022).

In this paper, we present the first application of this newly developed method to the two-dimensional velocity and magnetic vector fields, both confined to the (y, z) plane. The description of the method is provided in Sect. 2. In Sect. 3.1 we apply it to analytically prescribed solutions for the velocity and magnetic fields, whereas in Sect. 3.2 we apply it to a two-dimensional MHD simulation of the solar surface. A discussion on the errors in the inference of horizontal velocities in two dimensions is presented in Sect. 4. Finally, we present a summary of our findings in Sect. 5.

2. Numerical method

2.1. General remarks

Our ultimate goal is to infer the horizontal velocities v_x, v_y (i.e., those perpendicular to the line of sight direction). To this end, we used the induction equation under the assumption of ideal MHD (i.e., zero magnetic diffusion):

$$\begin{array}{c}\hfill \frac{\partial \mathbf{B}}{\partial t}=\mathrm{\nabla }\times (\mathbf{v}\times \mathbf{B}).\end{array}$$(1)

This equation describes how the magnetic field evolves over time in the presence of a velocity field and is satisfied at all times. Therefore, we can specify a particular instant in time, t₀:

$$\begin{array}{c}\hfill \frac{\partial \mathbf{B}}{\partial t}{|}_{t_0}=\dot{\mathbf{B}}(t_0)=\mathrm{\nabla }\times (\mathbf{v}(t_0)\times \mathbf{B}(t_0)).\end{array}$$(2)

As mentioned in Sect. 1, it is possible to observationally infer $\mathbf{B},\dot{\mathbf{B}}$ and v_z. Hence the solution of this equation will provide the horizontal components (i.e., perpendicular to the line of sight) of the velocity: v_x and v_y.

2.2. Two-dimensional induction equation in the (y, z) plane

Although the remarks in Sect. 2.1 apply to the general three-dimensional case, the goal of this paper is to first test our method in two dimensions, the (y, z) plane, where the z-direction corresponds to the observer’s line of sight. To this end, we set v_x = B_x = 0, and thus the horizontal velocity to be retrieved becomes v_y(y, z). Under these assumptions, the two nontrivial components of the induction equation can be written as

$$\begin{array}{cc}\hfill \dot{B_y}& =\frac{\partial }{\partial z}(v_y{B}_z-v_z{B}_y)\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \dot{B_z}& =\frac{\partial }{\partial y}(v_z{B}_y-v_y{B}_z),\hfill \end{array}$$(4)

where for simplicity t₀ is not indicated anymore. We further assumed that the Stokes inversion provides B_y, B_z, $\dot{B_y}$, $\dot{B_z}$, and v_z. With this, we obtained an overdetermined system with one unknown, v_y, and two equations for each point of the two-dimensional domain.

2.3. Discretization and matrix form

We used the finite differences method with a five-point centered stencil throughout the entire domain to solve the induction equation numerically. With this, the first derivative of a continuous function, f, along direction r at point i can be written as

$$\begin{array}{c}\hfill \frac{\partial f}{\partial r}{|}_i=\frac{f_{i-2}-8f_{i-1}+8f_{i+1}-f_{i+2}}{12h_r},\end{array}$$(5)

where h_r is the stencil length along direction r. Applying this stencil to Eqs. (3) and (4), we obtained

$$\begin{array}{cc}\hfill \dot{B_y}(i,j)+\frac{\partial }{\partial z}(v_z{B}_y){|}_{i,j}& =\frac{B_z(i,j)}{12h_z}{v}_y(i,j-2)\hfill \\ \hfill & -\frac{2B_z(i,j)}{3h_z}{v}_y(i,j-1)+\frac{\partial B_z}{\partial z}{|}_{i,j}{v}_y(i,j)\hfill \\ \hfill & +\frac{2B_z(i,j)}{3h_z}{v}_y(i,j+1)-\frac{B_z(i,j)}{12h_z}{v}_y(i,j+2)\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill \dot{B_z}(i,j)-\frac{\partial }{\partial y}(v_z{B}_y){|}_{i,j}& =-\frac{B_z(i,j)}{12h_y}{v}_y(i-2,j)\hfill \\ \hfill & +\frac{2B_z(i,j)}{3h_y}{v}_y(i-1,j)-\frac{\partial B_z}{\partial y}{|}_{i,j}{v}_y(i,j)\hfill \\ \hfill & -\frac{2B_z(i,j)}{3h_y}{v}_y(i+1,j)+\frac{B_z(i,j)}{12h_y}{v}_y(i+2,j),\hfill \end{array}$$(7)

where $\frac{\partial B_z}{\partial z}{|}_{i,j}$ and $\frac{\partial B_z}{\partial y}{|}_{i,j}$ are also evaluated according to Eq. (5). Equations (6) and (7) provide the basis to write an overdetermined linear system of equations where v_y(i, j) are the unknown quantities. The coefficients that multiply the unknown v_y in the equations above are renamed a_i, j or b_i, j, depending on whether they appear in the equations pertaining to $\dot{B_y}$ or $\dot{B_z}$, respectively. These coefficients are also given a superscript that refers to the grid point at which the horizontal component of velocity they multiply is evaluated. In addition, the independent coefficients (i.e., those that do not multiply v_y) are renamed $c_{i,j}^a$ and $c_{i,j}^b$. The superscript a or b indicates whether the independent coefficient appears in the equation for $\dot{B_y}$ or $\dot{B_z}$,. From Equation (6), we defined

$$\begin{array}{cc}\hfill a_{i,j}^{i,j-2}& =\frac{B_z(i,j)}{12h_z}\hfill \\ \hfill a_{i,j}^{i,j-1}& =-\frac{2B_z(i,j)}{3h_z}\hfill \\ \hfill a_{i,j}^{i,j}& =\frac{\partial B_z}{\partial z}{|}_{i,j}\hfill \\ \hfill a_{i,j}^{i,j+1}& =\frac{2B_z(i,j)}{3h_z}\hfill \\ \hfill a_{i,j}^{i,j+2}& =-\frac{B_z(i,j)}{12h_z}\hfill \\ \hfill c_{i,j}^a& =\dot{B_y}(i,j)+\frac{\partial }{\partial z}(v_z{B}_y){|}_{i,j}.\hfill \end{array}$$(8)

Likewise, for the b_i, j and c_i, j^b coefficients appearing in Equation (7),

$$\begin{array}{cc}\hfill b_{i,j}^{i-2,j}& =-\frac{B_z(i,j)}{12h_y}\hfill \\ \hfill b_{i,j}^{i-1,j}& =\frac{2B_z(i,j)}{3h_y}\hfill \\ \hfill b_{i,j}^{i,j}& =-\frac{\partial B_z}{\partial y}{|}_{i,j}\hfill \\ \hfill b_{i,j}^{i+1,j}& =-\frac{2B_z(i,j)}{3h_y}\hfill \\ \hfill b_{i,j}^{i+2,j}& =\frac{B_z(i,j)}{12h_y}\hfill \\ \hfill c_{i,j}^b& =\dot{B_z}(i,j)-\frac{\partial }{\partial y}(v_z{B}_y){|}_{i,j}.\hfill \end{array}$$(9)

With this, Eqs. (6) and (7) become

$$\begin{array}{cc}\hfill c_{i,j}^a& =a_{i,j}^{i,j-2}·v_y(i,j-2)+a_{i,j}^{i,j-1}·v_y(i,j-1)\hfill \\ \hfill & +a_{i,j}^{i,j}·v_y(i,j)+a_{i,j}^{i,j+1}·v_y(i,j+1)\hfill \\ \hfill & +a_{i,j}^{i,j+2}·v_y(i,j+2)\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill c_{i,j}^b& =b_{i,j}^{i-2,j}·v_y(i-2,j)+b_{i,j}^{i-1,j}·v_y(i-1,j)\hfill \\ \hfill & +b_{i,j}^{i,j}·v_y(i,j)+b_{i,j}^{i+1,j}·v_y(i+1,j)\hfill \\ \hfill & +b_{i,j}^{i+2,j}·v_y(i+2,j).\hfill \end{array}$$(11)

In matrix form, the system of equations can be written as $\widehat{\mathbf{A}}·\mathbf{x}=\mathbf{c}$, where $\widehat{\mathbf{A}}$ is constructed from the known coefficients a_i, j and b_i, j. This matrix has dimensions 2 ⋅ N ⋅ M × N ⋅ M for a domain of N × M grid points. x is the vector of unknowns with dimensions N ⋅ M × 1, while c corresponds to the vector of independent coefficients (c_i, j^a, c_i, j^b) with dimensions 2 ⋅ N ⋅ M × 1. In the matrix $\widehat{\mathbf{A}}$, we will have two rows (one for each available component of the induction equation used) associated with the same unknown v_y(i, j). The following piece of pseudo-code can be used to build the ${\widehat{A}}_{\mathit{kr}}$ elements of $\widehat{\mathbf{A}}$ as well as the c_k elements of the vector c:

1: r = 0

2: for k = 0, 2 ⋅ N ⋅ M − 1, 2 do

▹ Coordinates of the grid point (i, j) are given by

3: i = int(r/M)

4: j = mod(r, M)

 ▹  Elements referring to $\dot{B_y}$ (Eq. 3)

  $$\begin{array}{cc}\hfill A(k,i·M+j-2)& =a_{i,j}^{i,j-2}\hfill \\ \hfill A(k,i·M+j-1)& =a_{i,j}^{i,j-1}\hfill \\ \hfill A(k,i·M+j)& =a_{i,j}^{i,j}\hfill \\ \hfill A(k,i·M+j+1)& =a_{i,j}^{i,j+1}\hfill \\ \hfill A(k,i·M+j+2)& =a_{i,j}^{i,j+2}\hfill \\ \hfill c(k)& =c_{i,j}^a\hfill \end{array}$$(12)

 ▹   Elements referring to $\dot{B_z}$ (Eq. 4)

  $$\begin{array}{cc}\hfill A(k+1,(i-2)·M+j)& =b_{i,j}^{i-2,j}\hfill \\ \hfill A(k+1,(i-1)·M+j)& =b_{i,j}^{i-1,j}\hfill \\ \hfill A(k+1,i·M+j)& =b_{i,j}^{i,j}\hfill \\ \hfill A(k+1,(i+1)·M+j)& =b_{i,j}^{i+1,j}\hfill \\ \hfill A(k+1,(i+2)·M+j)& =b_{i,j}^{i+2,j}\hfill \\ \hfill c(k+1)& =c_{i,j}^b\hfill \end{array}$$(13)

5:  r = r + 1

6:end for

where the for-loop indexes refer to the starting value, end value, and step value. An example of the $\widehat{\mathbf{A}}$ matrix for a 10 × 10 domain is shown in Fig. 1. As explained above, for this particular example, $\widehat{\mathbf{A}}$ has a size of 200 × 100 (i.e., 200 equations for 100 unknowns). The color code indicates whether that particular point on the matrix comes from Eq. (3) (red) or Eq. (4) (blue).

Fig. 1. Structure of the matrix of the system $\widehat{\mathbf{A}}$ for a N × M = 10 × 10 domain defined by k rows and r columns. In red we show the position of the ${\widehat{A}}_{\mathit{kr}}$ coefficients related to Eq. (3) and in blue those related to Eq. (4).

2.4. Boundary conditions

To define the boundary conditions while maintaining the derivative scheme given by Eq. (5) throughout the entire domain, we introduced “ghost cells” around the boundaries of the domain (see Fig. 2). Although ghost cells are not part of the domain itself, we assumed a certain behavior of the unknown v_y function at those points. This behavior implies changes to the elements of the matrix $\widehat{\mathbf{A}}$ and to the components of the vector c related to the points close to the boundary. Since we used a five-point stencil scheme for the derivatives, the grid points whose coefficients will be modified are those with i = 0, 1, N − 2, N − 1 and/or j = 0, 1, M − 2, M − 1 (see Fig. 2), which refer to grid cells at the corners or sides of the domain, respectively.

Fig. 2. Example of a 2D grid for a 10 × 10 domain. The solution of the system is evaluated in the blue points. These points are called “inner cells”. The boundary condition must be specified in all red points. These are “ghost” or “boundary” cells. The white points in the corners are not used in the scheme considered here. In a different numerical scheme, they may play a role.

How the coefficients ${\widehat{A}}_{\mathit{rk}}$ of matrix $\widehat{\mathbf{A}}$ and the c_k components of vector c are modified close to the boundary depends on the actual boundary conditions employed. If we assume Dirichlet boundary conditions, we must set the value of v_y at the boundary. For example, v_y = f_D, where f_D denotes a fixed value and the subscript D indicates that Dirichlet boundary conditions are being applied. Consequently, the values of v_y in the ghost cells are no longer unknowns, and the terms of Eqs. (10) and (11) that involve v_y at the ghost cells become independent coefficients. For that reason, we get new c coefficients (referred to as c^a, ★ and c^b, ★) that can be written as a function of the old coefficients (c^a and c^b; given by Eqs. (8) and (9)) as

$$\begin{array}{cc}\hfill c_{0,0}^{a,★}& =c_{0,0}^a-a_{0,0}^{0,-2}·f_D(0,-2)-a_{0,0}^{0,-1}·f_D(0,-1)\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill c_{0,0}^{b,★}& =c_{0,0}^b-b_{0,0}^{-2,0}·f_D(-2,0)-b_{0,0}^{-1,0}·f_D(-1,0).\hfill \end{array}$$(15)

We note that in the case of Dirichlet boundary conditions, the ${\widehat{A}}_{\mathit{rk}}$ elements of matrix $\widehat{\mathbf{A}}$, built from a_i, j and b_i, j according to Eqs. (12) and (13), do not change.

On the other hand, we can assume symmetric boundary conditions by defining f₋₁ = f₀ and f₋₂ = f₁. These boundary conditions also ensure that, at some point in between the domain and the ghost zone, the derivative of v_y is zero (Neumann boundary condition). Unlike the case of the Dirichlet boundary conditions, where v_y = f_D was assumed to be known, Neumann boundary conditions link v_y in the ghost cells to the unknown v_y inside the domain. For this reason, close to the boundary we must modify the related ${\widehat{A}}_{\mathit{rk}}$ elements of matrix $\widehat{\mathbf{A}}$ given by Eqs. (8) and (9). The new coefficients, a^★ and b^★, derived from this assumption can be written as a function of the old a and b coefficients (given by Eqs. (12) and (13)). Here we provide examples for a corner point and a boundary line:

$$\begin{array}{c}\hfill \begin{array}{cc}{a}_{0,0}^{0,0,★}& =a_{0,0}^{0,0}+a_{0,0}^{0,-1}\\ a_{0,0}^{0,1,★}& =a_{0,0}^{0,1}+a_{0,0}^{0,-2}\\ b_{0,0}^{0,0,★}& =b_{0,0}^{0,0}+b_{0,0}^{-1,0}\\ b_{0,0}^{1,0,★}& =b_{0,0}^{1,0}+b_{0,0}^{-2,0}\end{array}\}\mathrm{if}(\mathrm{i},\mathrm{j})=(0,0)\end{array}$$(16)

$$\begin{array}{c}\hfill \begin{array}{cc}{b}_{0,j}^{0,j,★}& =b_{0,j}^{0,j}+b_{0,j}^{-1,j}\\ b_{0,j}^{1,j,★}& =b_{0,j}^{1,j}+b_{0,j}^{-2,j}\end{array}\}\mathrm{if}(\mathrm{i},\mathrm{j})=(0,\mathrm{j})\end{array}$$(17)

$$\begin{array}{c}\hfill \begin{array}{cc}{b}_{1,j}^{0,j,★}& =b_{1,j}^{0,j}+b_{1,j}^{-1,j}\end{array}\}\mathrm{if}(\mathrm{i},\mathrm{j})=(1,\mathrm{j}),\end{array}$$(18)

where 1 < j < n − 2 in Eqs. (17), (18). Finally, we note that in the case of Neumann boundary conditions, the c_k coefficients of c remain the same.

2.5. Solution of the linear system

Having defined the matrix of the system $\widehat{\mathbf{A}}$, the vector of independent coefficients c, and the boundary conditions, we could then solve the linear system of equations and infer the unknown horizontal velocity field v_y(y, z) that is encoded in the vector x. Owing to the fact that we have an overdetermined system with 2 ⋅ N ⋅ M equations and N ⋅ M unknowns, we solved the linear system in a least square sense: $\mathbf{x}={({\widehat{\mathbf{A}}}^{⊺}\widehat{\mathbf{A}})}^{-1}{\widehat{\mathbf{A}}}^{⊺}\mathbf{c}$, where the “^⊺” symbol indicates the transpose of the matrix. The solution of the linear system is found through the well-known LAPACK² package, in particular the subroutine dgesv that uses LU decomposition for a square matrix.

3. Comparison between analytical and numerical solutions

3.1. Two-dimensional analytical magnetic and velocity fields

We defined a two-dimensional domain in the (y, z) plane. This domain is discretized in N × M = 100 × 100 grid cells. Each grid cell is considered to be squared h_y = h_z = h with a length of h = 10⁶ cm. In this domain, we prescribe the following components for the magnetic and velocity fields at instant t₀:

$$\begin{array}{cc}\hfill B_y& =B_{1}exp\frac{\alpha }{\lambda }z\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill B_z& =B_{2}exp\frac{\alpha }{\lambda }y\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill v_y& =v_{1}cos\left\{\frac{2\pi (My-Nz)}{\mathit{NMh}}\right\}\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill v_z& =v_1\frac{M}{N}cos\left\{\frac{2\pi (My-Nz)}{\mathit{hNM}}\right\}+v_{2}cos\left\{\frac{2\pi y}{\mathit{hN}}\right\},\hfill \end{array}$$(22)

where B₁ = 350 Gauss, B₂ = 200 Gauss, α = −0.1, λ = 10h, v₁ = 10⁵ cm s⁻¹, and v₂ = v₁/2. We note that since B_y does not depend on y and B_z does not depend on z, the solenoidal condition for the magnetic field, ∇ ⋅ B = 0, is trivially satisfied. The detailed configuration of the magnetic field prescribed above is presented in Fig. 3.

Fig. 3. Two-dimensional analytical magnetic field. The color map shows the magnitude of the magnetic field ∥B∥. The green arrows show its direction. The color bar is saturated at 400 G to show the lower values of the magnetic field. The maximum value of the magnetic field at the selected region is ≈415 G.

Equations (19)–(22) were then applied to Equations (3)–(4) to find the analytical expressions for $\dot{B_y}$ and $\dot{B_z}$. In the following, we assumed $\dot{B_y}$, $\dot{B_z}$, B_z, B_y, and v_z to be known, leaving the horizontal component of the velocity as the only unknown, which we refer to as $\stackrel{\sim }{v_y}(y,z)$. To determine $\stackrel{\sim }{v_y}$, we employed vanishing Neumann boundary conditions instead of Dirichlet boundary conditions (see Sect. 2.4) because the latter require a knowledge of v_y at the boundaries, which in practical applications we do not have. To this end, we employed the expressions in Eqs. (16)–(18) extended to all the points close to the boundary (i.e., i = 0, 1, N − 2, N − 1 and j = 0, 1, M − 2, M − 1 for an N × M domain based on the indexing from Fig. 2).

With this, we applied the approach described in Sect. 2 to retrieve $\stackrel{\sim }{v_y}(y,z)$, which we then compared with the original v_y(y, z) given by Equation (21). Figure 4 shows the analytical velocity field v = (v_y, v_z) defined through Eq. (21) and Eq. (22) (panel a) and the inferred velocity field $\stackrel{\sim }{\mathbf{v}}=(\stackrel{\sim }{v_y},v_z)$ (panel b). Although at first glance it seems that both velocities are very similar, in the sense of direction and magnitude, one must bear in mind that these plots do not provide a complete picture because v_z is the same in both panels (i.e., it was assumed to be known).

Fig. 4. Velocity field maps in the plane (y, z). Panel (a) shows the velocity field defined analytically: v = (vy, vz) (see Eqs. (21) and (22)), whereas panel (b) provides the inferred velocity field: $\stackrel{\sim }{\mathbf{v}}=(\stackrel{\sim }{v_y},v_z)$. The modulus of the vectors is encoded in the color of the arrows as indicated in the color bar.

To further investigate the accuracy with which our method can determine the horizontal component of the velocity field in this analytical two-dimensional case, we show in Fig. 5 the relative error of the inferred velocity, defined as ${\epsilon }_{\stackrel{\sim }{v_y}}=|\frac{\stackrel{\sim }{v_y}-v_y}{v_y}|$ (left panel), and the angle formed by the analytical and the inferred velocity fields, defined as $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})=arccos(\mathbf{v}·\stackrel{\sim }{\mathbf{v}}/[\Vert \mathbf{v}\Vert \Vert \stackrel{\sim }{\mathbf{v}}\Vert ])$ (right panel).

Fig. 5. Errors of the inference of $\stackrel{\sim }{v_y}$ in the analytical 2D case. The left panel shows the relative error, ${\epsilon }_{\stackrel{\sim }{v_y}}$, of the inferred velocity. The right panel shows the angle between the analytical velocity (v) and the inferred one ($\stackrel{\sim }{\mathbf{v}}$). See text for more details. Both plots also indicate the mean and the median of these errors.

As can be seen, the average and median errors, ${\epsilon }_{\stackrel{\sim }{v_y}}$, over the entire domain are less than 1% and 1.5%, respectively. In addition, the angles between the original velocity field and the retrieved one, $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$, are also very small – about 0.5° on average. Over the majority of the domain, the actual retrieval is much better: ${\epsilon }_{\stackrel{\sim }{v_y}}\le 0.01\%$.

3.2. CO5BOLD numerical simulation in (y,z)

The analytical test in Sect. 3.1 features relatively simple magnetic and velocity fields (see e.g., Fig. 3). To test our method with a more realistic scenario, in this section we perform a similar test but using magnetic B(t₀) and velocity v(t₀) fields from a two-dimensional MHD numerical simulation of solar surface magneto-convection obtained from the CO5BOLD code (Freytag et al. 2012). The aforementioned simulation provides us with B_y(y, z),B_z(y, z),v_y(y, z), and v_z(y, z) at a fixed time step, t₀. In addition, the simulation has vanishing B_x and v_x. Details on how the 2D simulation was carried out are provided in Appendix A.

An important point to mention here pertains to the magnetic diffusivity, η. Although CO5BOLD solves the induction equation under ideal MHD, numerical diffusion is unavoidable, and therefore, if we were to take the magnetic and velocity fields at two different time steps, t₊ and t₋, separated by a small Δt = t₊ − t₋ and centered around t₀, the ideal induction equation would not be satisfied:

$$\begin{array}{c}\hfill \dot{\mathbf{B}}(t_0)={\displaystyle \underset{\mathrm{\Delta }t\to 0}{lim}}\frac{\mathbf{B}({\mathbf{t}}_{+})-\mathbf{B}({\mathbf{t}}_{-})}{2\mathrm{\Delta }t}\ne \mathrm{\nabla }\times [\mathbf{v}(t_0)\times \mathbf{B}(t_0)].\end{array}$$(23)

For this reason, we do not evaluate the time derivative at t₀ from two different snapshots (t₊ and t₋) of the MHD simulation. Instead we determine $\dot{\mathbf{B}}(t_0)$ directly from v(t₀) and B(t₀) using the right-hand side of Eq. (2). This ensures that the induction equation in ideal MHD is satisfied. This step can be justified because at this point we are only interested in studying the performance of the method described in Sect. 2 under ideal conditions. We defer the investigation of the effects of numerical or magnetic diffusivity, η_eff (see Appendix A), for a future study, which should also include the effect of having a finite Δt such that the limit in Eq. (23) is only approximately satisfied.

Next, we took a central region of N × M = 100 × 100 grid points from the full domain of the simulation and determined the horizontal component of the velocity v_y(y, z) from the known $\dot{\mathbf{B}}$, B, v_z and by employing vanishing Neumann boundary conditions. Figure 6 illustrates the magnetic field within this region, whereas Fig. 7 compares the original and the retrieved velocity fields. This last figure shows that the inferred velocity field is perceptually identical to the original from the MHD simulation.

Fig. 6. Two-dimensional CO5BOLD simulation magnetic field. The color map shows the magnitude of the magnetic field ∥B∥. The green arrows show its direction. The color bar is saturated at 1500 G to show the lower values of the magnetic field. The maximum value of the magnetic field within the selected region is ≈2245 G.

Fig. 7. Same as Fig. 4 but for the case of the two-dimensional MHD numerical simulation.

A more detailed study of the errors is presented in Fig. 8, where again the error in the inference of v_y, ${\epsilon }_{\stackrel{\sim }{v_y}}$ and the errors in the inference of the angle of the velocity field, $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$, are presented on the left and right panels, respectively. As can be seen, now the mean and median values of the errors are smaller than in the 2D analytical case (Sect. 3.1).

Fig. 8. Same as Fig. 5 but for the case of the two-dimensional MHD simulation.

4. Discussion

4.1. Errors

In the previous sections, we showed that the method developed in this work retrieves v_y(y, z) with a very high accuracy: about 1% mean error in the inference of v_y and less than 0.5° in the inference of the orientation of the total velocity field. We have referred to these two quantities as ${\epsilon }_{\stackrel{\sim }{v_y}}$ and $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$, respectively.

Interestingly, there are some specific regions where the values of ${\epsilon }_{\stackrel{\sim }{v_y}}$ and $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$ are much larger than the mean and median values. For instance, this is the case of the diagonal stripes in Fig. 5 (analytical case), where ${\epsilon }_{\stackrel{\sim }{v_y}}$ and $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$ reach values as high as 50% and 5°, respectively. Similar regions of large uncertainties also appear in Fig. 8 for the case of the 2D MHD simulation, albeit with smaller errors than in the analytical case.

To investigate the source of the large uncertainties in these regions, it is convenient to notice that, as per their definitions, ${\epsilon }_{\stackrel{\sim }{v_y}}$ will always be larger in those regions where v_y → 0, whereas $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$ will be larger in regions where ∥v∥→0. It is also important to note that ∥v∥→0 immediately implies v_y → 0, but not the other way around.

With this clarified, we are in a position to assert that in the analytical case (Sect. 3.1), ${\epsilon }_{\stackrel{\sim }{v_y}}$ and $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$ are the largest in regions where ∥v∥→0 (see Fig. 4). In the case of the two-dimensional MHD simulation (Sect. 3.2), the regions where the errors are largest do not coincide with locations where ∥v∥→0. We can ascertain this because in the MHD simulations, the minimum velocity is about 1 km s⁻¹ (see color bar in Fig. 7). In the regions where v_y → 0, and thus where the velocity is mostly aligned with the z-axis, ${\epsilon }_{\stackrel{\sim }{v_y}}$ instead increases. Also, the fact that the modulus of the velocity in the MHD simulations is never zero implies that the angle between the velocity vector in the MHD simulations and inferred velocity vector, $\angle (\mathbf{v},\stackrel{\sim }{\mathbf{v}})$, does not increase as much as in the analytical example presented in Sect. 3.1.

To showcase how well the horizontal component of the velocity v_y is inferred in a way that is independent from the magnitude of v_y or ∥v∥, we show in Fig. 9 scatter plots between the real v_y and the inferred values of the velocity $\stackrel{\sim }{v_y}$. One can see that in the analytical test (left; Sect. 3.1) and the test with MHD simulations (right; Sect. 3.2), the Pearson correlation coefficient (r) is almost 1.0. Furthermore, almost all points of the 2D domain fall almost exactly on top of the perfect retrieval line (i.e., slope one) indicated by the black line.

Fig. 9. Scatter plots of vy and $\stackrel{\sim }{v_y}$ for the analytical (left) and MHD numerical simulation (right) cases. The Pearson correlation coefficient is measured in both cases (r). The slope of one line (indicating a perfect retrieval) is displayed as a black line.

4.2. Known degeneracies

Due to the structure of the induction equation, there are theoretically two degeneracies that can affect the retrieval of the velocity vector. First, it is possible to add a gradient of a scalar function, ϕ, to the cross product v × B without affecting $\dot{\mathbf{B}}$. Second, only the velocity component perpendicular to the magnetic field contributes to its variation over time.

These two degeneracies suggest that the velocity field inferred using the method presented in this paper may not represent the real velocity of the plasma. This raises the question of how we can ensure that the retrieved velocity is the one we seek.

In two dimensions, neither degeneracy is very relevant. For instance, in the (y, z) plane, we define v and B both perpendicular to the x-axis. Then, v × B is parallel to the unit vector e_x. Moreover, if we consider a scalar function, ϕ, whose gradient can be added to the cross product, v × B, then by definition, ∇ϕ is also parallel to e_x. However, given that ϕ is only a function of (y, z), it trivially follows that $\mathrm{\nabla }\varphi =\frac{\partial \varphi (y,z)}{\partial x}{\mathbf{e}}_x=0$.

In three dimensions, the two aforementioned degeneracies can potentially play a much more important role. We anticipate being able to retrieve the real plasma velocity field in 3D for two reasons. First, v_z is considered to be known (i.e., from the Stokes inversion that measures it from the Doppler effect), and this component of the velocity already includes contributions from the parallel and perpendicular components of the velocity with respect to B. For this reason, the knowledge of v_z breaks the degeneracy, and it allows us to retrieve the full velocity field through the induction equation. Second, our method can be generalized to include other MHD constraints, such as mass conservation. These additional constraints could help us obtain the plasma velocity field independently of the term ∇ϕ that can be added to v × B. We will address these topics in more detail in future work.

5. Conclusions and future work

We have developed a new method to infer the velocities perpendicular to the line of sight in the solar atmosphere. This method employs the full induction equation in ideal MHD and assumes that the line-of-sight component of the velocity, v_z; the magnetic field, B; and its time derivative, $\dot{\mathbf{B}}$, are known as a function of space and time. Nowadays, this can be achieved thanks to novel Stokes inversion techniques applied to data from modern spectropolarimeters (Kaithakkal et al. 2020; Liu et al. 2025).

As a consequence of using the three components of the induction equation, the method also allows retrieval of the full z dependence of the horizontal velocity. This is not possible when employing only the vertical (z direction) component of the induction equation (see, e.g., Longcope 2004).

As a first step, we tested our newly developed method in two-dimensions by first using analytically prescribed functions for the velocity and magnetic fields and secondly by using velocity and magnetic fields from MHD simulations of the solar surface magnetoconvection. In these 2D cases, only two of the components of the induction equation were used. Our results indicate that in an ideal scenario, it is possible to retrieve the horizontal component of the velocity with an average error of about 1%. In regions where the horizontal component of the velocity is very small, uncertainties in the direction and relative error in the magnitude become significantly larger. Overall, Pearson’s correlation coefficient is very close to unity (r ≃ 1) in both of the studied cases. This allows us to conclude that our method can retrieve the velocity correctly.

An important feature of the method presented in this paper is that it can be easily extended to three dimensions by including the third component of the induction equation. Although the discretization of the equations and the construction of the matrix $\widehat{\mathbf{A}}$ follow the same procedure described in this paper (see Eqs. (12)–(13)), the matrix dimensions increase substantially as the domain expands from two to three dimensions. Specifically, the size of the matrix, $\widehat{\mathbf{A}}$, changes from 2NM × NM to 3NML × 2NML because the system of equations includes an additional dimension as well as an additional equation with two unknowns (v_x, v_y) at each point of the domain. Inverting a matrix of this size requires a considerable amount of computational time, but we plan on undertaking this step next. Apart from extending our approach to three dimensions, in the future we hope to investigate the accuracy with which we can infer the horizontal components of the velocity under more realistic scenarios. Namely, we plan to examine the presence of random or systematic errors in the input parameters (v_z, B, and $\dot{\mathbf{B}}$) given by the Stokes inversion, the effects of having a finite time sampling when determining $\dot{\mathbf{B}}$ observationally, and the existence of a finite magnetic diffusion. All of these additional sources of uncertainty will likely overshadow the 1% uncertainty of the idealized scenarios investigated here.

Acknowledgments

The work presented in this paper has been funded by a grant from the Deutsche Forschung Gemeinschaft (DFG): project 538773352. This research has made use of NASA’s Astrophysics Data System. The authors are grateful to Petri Käpylä, Maria Kazachenko, Brian Welsch, and Oskar Steiner for useful discussions on the topic.

References

Appendix A: Two-dimensional CO5BOLD simulations

We used the CO5BOLD code (Freytag et al. 2012) in the box-in-a-star setup to create a realistic 2D atmosphere. The code solves the equations of ideal MHD for a fully compressible gas in a Cartesian box using a finite-volume formulation with an approximate Riemann type solver (HLLMHD; Harten et al. 1983; Schaffenberger et al. 2005; Steiner et al. 2013). The constrained transport method used in CO5BOLD ensures the magnetic field divergence free condition (∇ ⋅ B = 0) to machine precision. Although, CO5BOLD solves an ideal MHD induction equation, the effective magnetic diffusivity due to the discretization is found to be around η_eff ≈ 3.0 × 10¹⁰ cm² s⁻¹ for a computational cell width of 10 km and at the bottom of the photosphere (Riva & Steiner 2022). For our model, we used a realistic solar-like equation of state taken from the CIFIST project (Ludwig et al. 2009). The radiative transfer proceeds via a short characteristics scheme as described in Steffen (2017), with gray opacities from the MARCS model atmosphere package (Gustafsson et al. 2008), provided in tabulated form as a function of gas pressure and temperature.

For this work, we selected a snapshot from a 3D magnetoconvection simulation computed using CO5BOLD, which was originally obtained by embedding a homogeneous, vertical field of 50 G magnetic flux density in a relaxed model of convection (model d3gt57g44v50fc from Calvo 2018). We extracted a 2D slice in the (y, z) plane from this snapshot at a location in the x-direction that contains a few kilogauss flux concentrations. The selected computational domain spans 9.6 Mm x 2.8 Mm in the (y, z) plane, with a cell size of 10 km in both spatial directions, of which a smaller domain is selected for our study (see Fig. 6). In CO5BOLD, all the quantities are defined at cell centers, except for the magnetic field components, which are defined at the centers of their corresponding faces. A constant external gravity field with g = 275 m s⁻² acts along the z-direction. The top boundary is open for fluid flow and outward radiation, with the density decreasing exponentially in the boundary cells outside the domain. The vertical component of the magnetic field (B_z) is constant across the top boundary, while the transverse component (B_y) drops to zero at the same location.

The bottom boundary conditions for the magnetic fields are the same as for the top boundary. For the side boundaries, we used periodic boundary conditions, except for the magnetic field in the x coordinate direction. For this direction, the x component of the magnetic field is set to zero, and a constant extrapolation applies to the other two components. We note that even though initially B_x and v_x are different from zero, the simulation is advanced for an additional 1000 seconds of solar time, after which B_x and v_x become negligible.

The bottom boundary for the thermodynamic quantities is set up in such a way that the in-flowing material carries a constant specific entropy of 1.775 × 10⁹ erg g⁻¹ K⁻¹ resulting in a radiative flux corresponding to an effective temperature (T_eff) of ∼5770 K.

All Figures

	Fig. 1. Structure of the matrix of the system $\widehat{\mathbf{A}}$ for a N × M = 10 × 10 domain defined by k rows and r columns. In red we show the position of the ${\widehat{A}}_{\mathit{kr}}$ coefficients related to Eq. (3) and in blue those related to Eq. (4).
In the text

	Fig. 2. Example of a 2D grid for a 10 × 10 domain. The solution of the system is evaluated in the blue points. These points are called “inner cells”. The boundary condition must be specified in all red points. These are “ghost” or “boundary” cells. The white points in the corners are not used in the scheme considered here. In a different numerical scheme, they may play a role.
In the text

	Fig. 3. Two-dimensional analytical magnetic field. The color map shows the magnitude of the magnetic field ∥B∥. The green arrows show its direction. The color bar is saturated at 400 G to show the lower values of the magnetic field. The maximum value of the magnetic field at the selected region is ≈415 G.
In the text

	Fig. 4. Velocity field maps in the plane (y, z). Panel (a) shows the velocity field defined analytically: v = (vy, vz) (see Eqs. (21) and (22)), whereas panel (b) provides the inferred velocity field: $\stackrel{\sim }{\mathbf{v}}=(\stackrel{\sim }{v_y},v_z)$. The modulus of the vectors is encoded in the color of the arrows as indicated in the color bar.
In the text

Fig. 5. Errors of the inference of $\stackrel{\sim }{v_y}$ in the analytical 2D case. The left panel shows the relative error, ${\epsilon }_{\stackrel{\sim }{v_y}}$, of the inferred velocity. The right panel shows the angle between the analytical velocity (v) and the inferred one ($\stackrel{\sim }{\mathbf{v}}$). See text for more details. Both plots also indicate the mean and the median of these errors.

In the text

	Fig. 6. Two-dimensional CO5BOLD simulation magnetic field. The color map shows the magnitude of the magnetic field ∥B∥. The green arrows show its direction. The color bar is saturated at 1500 G to show the lower values of the magnetic field. The maximum value of the magnetic field within the selected region is ≈2245 G.
In the text

	Fig. 7. Same as Fig. 4 but for the case of the two-dimensional MHD numerical simulation.
In the text

	Fig. 8. Same as Fig. 5 but for the case of the two-dimensional MHD simulation.
In the text

	Fig. 9. Scatter plots of vy and $\stackrel{\sim }{v_y}$ for the analytical (left) and MHD numerical simulation (right) cases. The Pearson correlation coefficient is measured in both cases (r). The slope of one line (indicating a perfect retrieval) is displayed as a black line.
In the text