Magnetic reconnection, a fundamental plasma process, is pivotal in understanding energy conversion and particle acceleration in astrophysical systems. While extensively studied in two-dimensional (2-D) configurations, the dynamics of reconnection in three-dimensional (3-D) systems remains under-explored. In this work, we extend the classical tearing mode instability to three dimensions by introducing a modulation along the otherwise uniform direction in a 2-D equilibrium, given by g(y), mimicking a flux-tube-like configuration. We perform linear stability analysis (both analytically and numerically) and direct numerical simulations to investigate the effects of three-dimensionality. Remarkably, we find that a tearing-like instability arises in three dimensions as well, even without the presence of guide fields. Further, our findings reveal that the 3-D tearing instability exhibits reduced growth rates compared with two dimensions by a factor of ∫ g(y)^(1/2) dy / ∫ dy, with the dispersion relation maintaining similar scaling characteristics. We show that the modulation introduces spatially varying resistive layer properties, which influence the reconnection dynamics. Magnetic reconnection is a fundamental plasma process in the contexts of astrophysical, space and fusion plasmas (Zweibel & Yamada 2009). This process occurs when oppositely directed magnetic field lines interact in the regions of high current density to break and reconnect, altering the magnetic topology and allowing for the rapid conversion of magnetic energy into kinetic and thermal energy. Magnetic reconnection can drive explosive phenomena such as solar flares, coronal mass ejections and geomagnetic storms in the Earth’s magnetosphere (Shibata & Magara 2011; Burch & Phan 2016; Ruan, Xia & Keppens 2020). It has played a pivotal role in regulating the dynamics of high-energy astrophysical environments, from the solar corona and interstellar medium to distant pulsar magnetospheres and black hole accretion disks (Cerutti, Philippov & Spitkovsky 2016; Ripperda, Bacchini & Philippov 2020; Fielding, Ripperda & Philippov 2023; Zhang, Pree & Bellan 2023). In particular, reconnection is invoked to understand the particle acceleration and non-thermal emission in many systems (Cerutti, Uzdensky & Begelman 2012; Sironi & Spitkovsky 2014; Guo et al. 2015; Brunetti & Lazarian 2016; Werner et al. 2016; Ghosh & Bhat, 2025). Uncovering the mechanisms of reconnection is essential for explaining these high-energy events and for improving our understanding of plasma behavior across the universe. Reconnection has been largely studied in two dimensions. It is thought to manifest in different ways: (i) spontaneously as the tearing mode instability, (ii) in steady state known as the Sweet–Parker (SP) model and (iii) as turbulent reconnection. Tearing instability was first studied in the context of magnetic confinement in laboratory plasmas (Furth, Killeen & Rosenbluth 1963; Coppi et al. 1976). Since both tearing modes and the SP model lead to dimensionless reconnection rates that depend on a negative power of the Lundquist number (S), they could not explain observed time scales pertaining to solar flares. However, more recently, the discovery of the plasmoid instability has been considered to have solved the time scale problem (Loureiro, Schekochihin & Cowley 2007; Bhattacharjee et al. 2009; Cassak, Shay & Drake 2009; Uzdensky, Loureiro & Schekochihin 2010; Pucci & Velli 2013; Comisso et al. 2016). The plasmoid instability is fundamentally a re-emergence of the tearing instability in high-S regimes and arises asymptotically beyond S ~ 10^4, leading to bursty reconnection and the formation of small secondary magnetic islands called plasmoids (Samtaney et al. 2009; Landi et al. 2015). Importantly, in its nonlinear regime, the plasmoid instability yields reconnection rates that are effectively independent of S (Bhattacharjee et al. 2009). Complementary to this, the ideal tearing scenario proposed by Pucci & Velli (2013) predicts that thin current sheets can become unstable on ideal (Alfvénic) time scales when their aspect ratio scales as S^(-1/3), leading to reconnection rates that are independent of resistivity. This prediction has been confirmed and extended through simulations and analyses in various geometries and conditions (Landi et al. 2015; Del Zanna et al. 2016; Papini, Landi & Zanna 2019). Alternatively, there has been work that shows that three-dimensional (3-D) turbulence can also help in leading to fast reconnection (Lazarian et al. 2020). This turbulent model intrinsically uses the SP model for understanding the local reconnections and thus it is unclear if it needs some modification given that the SP model is ruled out at values of S higher than ~10^4. In this work, rather than focusing on plasmoid instability or turbulent reconnection models, we aim to extend the study of the tearing mode instability into a fully 3-D context. This provides a complementary and necessary perspective for understanding reconnection in more realistic, inherently 3-D systems. This has been approached in several ways, with two of the simplest being : (i) extending a 2-D initial equilibrium into the third dimension and introducing 3-D perturbations, (ii) incorporating a uniform guide field along the third dimension. Approach (i) has been explored to demonstrate the occurrence of the kink instability, where the equilibrium current sheet buckles in response to 3-D perturbations (Landi et al. 2008; Oishi et al. 2015). This buckling leads to nonlinear reconnection processes that can be faster than their 2-D counterparts. Reconnection set-ups are sensitive, even during the linear phase, to parameters such as whether the initial equilibrium is modeled using pressure balance or force-free fields, as well as the inclusion of a guide field (Landi et al. 2008). In the study by Onofri et al. (2004), which falls under category (ii), the presence of a guide field was found to stabilize the 3-D instabilities observed by Dahlburg, Antiochos & Zang (1992), resulting in behavior that is closer to a quasi-2-D dynamics. They also observed that, in the nonlinear regime, island coalescence leads to faster reconnection compared with the linear regime. Faster reconnection in the nonlinear phase appears to be a recurring finding across 3-D tearing set-ups. For instance, Wang, Yokoyama & Isobe (2015) showed that introducing random perturbations, rather than specific mode perturbations, results in the formation of multiple tearing layers. These layers interact, leading to faster reconnection. Ultimately, many of these studies report that extending reconnection set-ups into three dimensions facilitates the generation of turbulence, which plays a crucial role in enhancing the reconnection dynamics. Another approach, taken primarily by solar physicists, is to investigate 3-D field configurations that either have null points (where the magnetic field strength vanishes) or possess layers conducive to reconnection (Parnell et al. 2010). Studies of magnetic null points have identified distinct topological features, such as spine lines and fan surfaces (Priest & Démoulin 1995). These structures enable reconnection by directing magnetic flux along separatrix surfaces, which divide regions of differing magnetic connectivity (Wyper & Pontin 2014). However, many reconnection events, particularly in solar and astrophysical plasmas, occur in regions lacking null points (Démoulin et al. 1997) In these cases, quasi-separatrix layers (QSLs) are thought to play a crucial role. The QSLs are regions where magnetic field lines experience rapid connectivity changes, even without intersecting at null points. It is proposed that reconnection within QSLs can occur across a distributed region in the presence of intense current layers and, under certain conditions, exhibit bursty, explosive behavior similar to null-point reconnection (Aulanier et al. 2006; Baker et al. 2009; Kumar et al. 2021; Mondal et al. 2023). In this work, we explore the reconnection between anti-parallel, flux-tube-like fields. This draws inspiration from the configuration used to examine vortex tube reconnection in Melander & Hussain (1989), where the tubes, with cylindrical symmetry, are characterized by an axial component of the field with a radial dependence. We adopt a simplified version of this set-up, which can be considered a modulation of the classic Harris sheet in the third dimension, described by B_z(x) = tanh(x) sech^2(x) sech^2(y). Flux tube reconnection has been studied commonly with an intent to mainly explore the effect of twists and writhes in the field on the ensuing reconnection (Dahlburg & Antiochos 1997; Linton, Dahlburg & Antiochos 2001; Wilmot-Smith & De Moortel 2007). These configurations involving significant helicity are often referred to as flux ropes and exhibit a complex 3-D dynamics. In contrast, flux tubes without helicity are simpler, providing an idealized framework for examining fundamental aspects of magnetic reconnection. However, non-helical flux-tube reconnection has received less attention in the literature compared with flux ropes (Linton & Priest 2003). In both cases – whether with or without helicity – previous studies have primarily considered flux tube interactions with either finite inclination angles or perpendicular orientations relative to each other. These interactions typically involve flattening of the tubes and formation of topologically complex structures during reconnection. In contrast, our study focuses on the simpler and less explored case of zero inclination angle, providing new insights into this idealized configuration. The further structure of the paper is as follows. In § 2, we describe the linear stability analysis using analytical and numerical approaches to understand the effects of modulation along the third dimension. Section 3 outlines the set-up of the direct numerical simulations used to test the theory and obtain the tearing mode growth rates, and § 4 presents the results from the simulations comparing them with 2-D cases and discussing the impact of 3-D effects. Finally, § 5 summarizes our findings and suggests areas for future research. To investigate 3-D effects in the tearing instability, we consider a 3-D base state that consists of a modulation of the standard 2-D equilibria, given by a magnetic field configuration of the form B_0 = B_0 f(x) z_hat, in the third direction. The 3-D initial configurations are obtained by modulating the corresponding 2-D configurations along the third direction by a function g(y), such that the 3-D base state is then given by B_0 = B_0 f(x) g(y) z_hat. To have a configuration of reversing magnetic fields along x, f(x) is chosen to be an odd function. Moreover, as we shall show in § 3, a smoothly varying g(y) with even parity such that g(y) -> 0 as |y| -> infinity produces a tube-like configuration of magnetic fields. We choose such a g(y). Since the magnetic field in such a configuration varies only across (and not along) itself, the associated magnetic tension, B · ∇B = 0. The initial gas pressure distribution can then be chosen to cancel the magnetic pressure, thus ensuring that this initial configuration is in equilibrium. Here, we have ignored the dissipation of the magnetic field, as is usually done in the tearing mode derivation under the assumption that the diffusive time scales are much larger than the tearing time scales. We consider the full incompressible and inviscid magnetohydrodynamic (MHD) equations: (2.1) ∇ · u = 0, ∇ · B = 0 (2.2) ∂(∇ × u)/∂t = - u · ∇(∇ × u) + ∇ × [(∇ × B) × B] (2.3) ∂B/∂t = ∇ × (u × B) + η ∇^2 B which comprise the solenoidality condition for the velocity, u, and the magnetic field, B, and the evolution equations for the vorticity (from the momentum equation) and the magnetic field (the induction equation), respectively. Here, η is the magnetic diffusivity, which quantifies the rate at which magnetic field lines diffuse through the plasma due to finite electrical conductivity. We work in Alfvén units wherein 4πρ = 1. The full MHD equations are linearized about the base state B_0 f(x) g(y) z_hat in terms of the perturbed magnetic fields, b, and the perturbed velocity, u. Assuming perturbed quantities of the form ψ = ψ(x,y) exp(i k z - i ω t), the linearized vorticity and induction equations yield: (2.7) -iω (k^2 u_y - ∂_x ∂_y u_x - ∂_y^2 u_y) = i k B_0 [f(x) g''(y) b_y + f'(x) g'(y) b_x + f'(x) g(y) ∂_y b_x - f(x) g'(y) ∂_x b_x - f(x) g(y) ∂_y ∂_x b_x - f(x) g(y) ∂_y^2 b_y + f(x) g(y) k^2 b_y] (2.8) -iω (∂_x^2 u_x + ∂_x ∂_y u_y - k^2 u_x) = i k B_0 [-f'(x) g'(y) b_y - f(x) g'(y) ∂_x b_y - f''(x) g(y) b_x + f'(x) g(y) ∂_y b_y - f(x) g(y) k^2 b_x + f(x) g(y) ∂_x^2 b_x + f(x) g(y) ∂_x ∂_y b_y] (2.9) -iω b_x = i k B_0 f(x) g(y) u_x + η (∂_x^2 b_x + ∂_y^2 b_x - k^2 b_x) (2.10) -iω b_y = i k B_0 f(x) g(y) u_y + η (∂_x^2 b_y + ∂_y^2 b_y - k^2 b_y) We adopt the approach in Goldston & Rutherford (1995) for further analysis. The current sheet is divided into three regions: inner, outer, and overlap regions. The instability parameter is defined as: (2.11) Δ' = [∂ ln(b_x)/∂x] from 0^- to 0^+, where b_x = b_x(x,y). In the inner region, applying standard considerations ∂_x^2 >> k^2, f(x) ≈ x, and ∂_x^2 >> ∂_y^2: (2.12) ∂_x^2 b_x = (-1/η) (iω b_x + i k B_0 x g(y) u_x) (2.13) (-ω/k) (∂_x^2 u_x + ∂_x ∂_y u_y) = B_0 [-x g(y) ((iω b_x + i k B_0 x g(y) u_x)/η) - g'(y) b_y - x g'(y) ∂_x b_y + g(y) ∂_y b_y + x g(y) ∂_x ∂_y b_y] Balancing terms yields the characteristic resistive layer width: (2.14) x ~ δ = (η γ)^(1/4) / (B_0 k g(y))^(1/2), where γ = -iω. The velocity scale is: (2.15) u_x_tilde ~ i γ b_x_tilde / (k B_0 g(y) δ). Integrating over the inner layer gives: (2.16) [∂_x b_x] = (1/η) ∫ (γ b_x_tilde^2 - i k B_0 x g(y) u_x) dx (2.17) (1/b_x_tilde) [∂_x b_x] = (γ δ(y)/η) ∫ (1 + X V) dX = (γ δ(y)/η) I(y) (2.18) Δ' = (γ/η) ((γ η)^(1/4) / (k B_0 g(y))^(1/2)) I(y) Integrating both sides over y to eliminate the y-dependence for a global growth rate γ gives: (2.19) Δ' ∫ g(y)^(1/2) dy = (γ^(5/4) / (η^(3/4) (k B_0)^(1/2))) ∫ I(y) dy Rearranging gives the 3-D dispersion relation: (2.20) γ = (Δ'^(4/5) η^(3/5) (k B_0)^(2/5) ∫ g(y)^(1/2) dy) / (∫ I(y) dy) This shows that the growth rate in 3D is reduced compared to 2D by a factor of ∫ g(y)^(1/2) dy / ∫ dy, while maintaining similar scaling characteristics. The linearized equations (2.7)–(2.10) can be formulated as a generalized eigenvalue problem: (2.21) γ M v = L v where M and L are linear operators acting on the state vector v = [u_x, u_y, b_x, b_y], and γ = -iω is the eigenvalue. We choose a modified Harris sheet as the unmodulated base state: f(x) = 2.6 tanh(x) sech^2(x), where the prefactor 2.6 ensures a peak magnetic field strength of 1 (v_A = 1). With modulation along y, the 3-D base state is: (2.22) B = f(x) g(y) z_hat = 2.6 tanh(x) sech^2(x) sech^2(y/λ) z_hat where λ controls the modulation width and the shear length is set to a = 1. Numerical solution via pseudo-spectral methods and Fourier decomposition (using PyTorch and SciPy) yields the 3D eigenfunctions and dispersion curves. The 3D eigenfunctions display converging velocity fields at the origin, a double-humped profile for b_x, and a quadrupolar structure for b_y. The dispersion curves confirm that 3D growth rates are uniformly attenuated by the modulation while preserving the overall dispersion shape. We perform direct numerical simulations in two and three dimensions to obtain the corresponding dispersion relations using a pseudo-spectral code written with the Dedalus framework. Initial equilibrium is maintained by balancing magnetic pressure with gas pressure (p_gas + p_mag = const). The Lundquist number is defined as S = v_A a / η and Alfvénic time scale is τ_A = a / v_A. To measure mode-by-mode growth rates, single-wavenumber perturbations b_x = sin(k* z) are introduced in a domain of size L_x = 4π, L_y = 4π, L_z = 2π/k* with grid resolutions N_x = 128, N_y = 128, N_z = 64. To prevent background equilibrium diffusion (-η ∇^2 B_0) from decaying the base state on timescales comparable to the tearing growth at moderate Lundquist numbers, a compensating source term is added to the induction equation. Direct numerical simulations demonstrate the full evolution of the 3-D tearing instability. The anti-parallel flux tubes converge, break, and reconnect to form new flux-tube structures. The growth rate is calculated by tracking the spectral energy E_k*(b) = ∑ |b_hat(k_x, k_y, k*)|^2 over time in the linear growth phase. Measurements across different values of λ and k show excellent agreement with both the numerical eigenvalue solver and the analytical dispersion relation (2.20). When rescaled by the predicted prefactor ∫ g(y)^(1/2) dy / ∫ dy, the 3D dispersion curves collapse exactly onto the 2D curve. The scaling of the maximum growth rate with the Lundquist number follows the classic FKR/tearing scaling, γ_max ~ S^(-1/2), across simulations with S = 80, 160, and 320. Tracking the perturbation growth across individual y-slices reveals identical growth rates, demonstrating cross-coupling across the 3D domain rather than uncoupled 2D slice dynamics. The 3D magnetic islands exhibit a bent structure in the y-direction toward the enhanced current density region around the central reconnection site. We investigated a 3-D extension of the tearing instability with modulated flux-tube-like magnetic fields without guide fields. We found that tearing-like modes arise naturally in this configuration, amenable to tractable linear theory. The linear growth rate is reduced by the factor ∫ g(y)^(1/2) dy / ∫ dy due to non-uniform resistive layer thickness along the third dimension. Reconnection occurs plane by plane with cutting and rejoining of field lines while maintaining a coherent 3-D structure. Crucially, the maximum growth rate scaling γ_max ~ S^(-1/2) is preserved in 3D even without a guide field. Future avenues of investigation include the manifestation of the plasmoid instability in such configurations, the influence of magnetic helicity, turbulent background interactions, and quasi-separatrix layer dynamics. A.1. Assumption of ∂_x^2 >> ∂_y^2: Verified from the calculated eigenfunctions; the ratio |∂_x^2 b / ∂_y^2 b| is significantly greater than 1 in the region of high magnetic shear. A.2. Details of the Eigenvalue Solver: Linear operators M and L are expressed as sparse matrix operators of size (4 N_x N_y × 4 N_x N_y) in Fourier basis, and solved for eigenvalues and eigenvectors using PyTorch GPU acceleration. The codebase is publicly available as SPEC-Tear. A.3. Low-pass filtering and periodicity: High-frequency Fourier modes of the non-periodic analytic equilibrium profiles are truncated above a threshold to ensure strict periodicity in the spectral solver. A.4. Magnetic island structure and null points: Visualizations of magnetic field lines confirm non-zero b_y, producing bent island structures with symmetric X-point reconnection sites across the y=0 plane.