Global Chemical Transport on Hot Jupiters: Insights from 2D VULCAN photochemical model

The atmospheric dynamics of tidally-locked hot Jupiters is dominated by the equatorial winds. Understanding the interaction between global circulation and chemistry is crucial in atmospheric studies and interpreting observations. Two-dimensional (2D) photochemical transport models shed light on how the atmospheric composition depends on circulation. In this paper, we introduce the 2D photochemical transport model, VULCAN 2D, which improves on the pseudo-2D approaches by allowing for non-uniform zonal winds. We extensively validate our VULCAN 2D with analytical solutions and benchmark comparisons. Applications to HD 189733 b and HD 209458 b reveal distinct characteristics in horizontal transport-dominated and vertical mixing-dominated regimes. Motivated by the inferred carbon-rich atmosphere by Giacobbe et al. (2021), we find that HD 209458 b with super-solar carbon-to-oxygen ratio (C / O) exhibits pronounced C 2 H 4 absorption on the morning limb but not on the evening limb, owing to horizontal transport from the nightside. We discuss when a pseudo-2D approach is a valid assumption and its inherent limitations. Finally, we demonstrate the e ff ect of horizontal transport in transmission observations and its impact on the morning-evening limb asymmetry with synthetic spectra, highlighting the need to consider global transport when interpreting exoplanet atmospheres.


INTRODUCTION
We are entering the era of detailed exoplanet atmosphere characterization.The atmospheric characterization has come a long way since the first transit observation with the Hubble Space Telescope (HST) (Charbonneau et al. 2002).Recent JWST spectral measurements revealed unprecedented details of gas giants (e.g., JWST Transiting Exoplanet Community Early Release Science Team et al. 2023;Rustamkulov et al. 2023;Bell et al. 2023).Transmission spectroscopy is central for probing atmospheric composition.To improve our ability to interpret observations through theoretical models, it is important to consider how the temperature and composition at the terminators and various regions are shaped by global circulation.This aspect cannot be addressed by traditional 1D models that omit thermal and chemical variations and transport processes in the horizontal direction.
Since transmission observations probe terminator regions that are composed of opposite sides of the planet, recent works have highlighted the need to account for spatial inhomogeneities.Espinoza & Jones (2021) and Grant & Wakeford (2023) introduced methods to separate the atmospheric components from the morning (leading) and evening (trailing) limbs.The advancement of high-resolution spectroscopy also expands our capacity to resolve the climate (e.g., Louden & Wheatley 2015) and chemical variations (e.g., Ehrenreich et al. 2020;Wardenier et al. 2021) across the planet.Atmospheric retrieval frameworks are beginning to extend to 2D and 3D to prevent biases from fitting the data with 1D profiles (MacDonald et al. 2020;Chubb & Min 2022;Pluriel, William et al. 2022).As we gain more complete observations from transit, eclipse, and phase curves, we will soon have the ability to probe the composition distribu-tions across the planet.It is essential to have models resolving global variations to support the progress of observations.
The atmospheric composition of gas giants is generally governed by thermochemistry, photochemistry, and mixing processes.The merit of 1D models is they can incorporate all the above processes as needed.However, the 1D structure intrinsically neglects the variations in 3D, and more importantly, how the global circulation impacts the local properties (Mendonc ¸a et al. 2018a;Drummond et al. 2020;Zamyatina et al. 2023).3D general circulation models (GCMs) can capture the intricate dynamic interactions across the planet and have become commonly applied to interpret observations (e.g., Beltz et al. 2021;Kempton et al. 2023).However, simplifications of the physical and chemical processes are usually required due to computation limitations.In addition to reducing the full Navier-Stokes dynamics into primitive equations (Mayne et al. 2019), 3D GCMs often make further simplifications to radiative transfer and chemical processes.The chemical processes in particular are often severely simplified in most applications, using the assumption of thermochemical equilibrium, which is valid only at high temperatures and pressures.Modeling efforts have been put into implementing reduced chemical schemes in a GCM to account for transport-induced disequilibrium (Drummond et al. 2020;Tsai et al. 2022;Lee et al. 2023).While Earth climate models have provided insights into the oxygen response of Earth-like atmospheres (Chen et al. 2021;Cooke et al. 2022;Braam et al. 2023;Deitrick & Goldblatt 2023;Ji et al. 2023), to date, the incorporation of photochemistry into 3D models for diverse non-Earth-like atmospheres has not been realized yet.
For most tidally-locked giant exoplanets, the stationary day-night heating drives large-scale equatorial waves.These waves interact with background flows and transport momentum from higher latitudes toward the equator, leading to an equatorial superrotating jet (Showman & Polvani 2011;Tsai et al. 2014;Showman et al. 2020).The strong jet plays a central role in global circulation, facilitating the transport of heat, aerosols, and gaseous species (Komacek & Showman 2016;Hammond & Lewis 2021;Steinrueck et al. 2021;Drummond et al. 2020).For these planets, "pseudo-2D" photochemical models (Agúndez et al. 2014;Venot et al. 2020;Baeyens et al. 2022;Moses et al. 2022) that employ a rotating 1D column to mimic a uniform jet in a Lagrangian frame have emerged as useful complementary tools.Their relatively fast computation enables them to incorporate the same detailed mechanisms as 1D models, but the simulated circulation is limited to uniform winds1 .On the other hand, 2D models implemented with horizontal diffusion have been used to model the meridional plane for solar system gas giants (Liang et al. 2005;Zhang et al. 2013;Hue et al. 2018) but have not yet been applied to exoplanets.In this work, we present the 2D (horizontal and vertical) photochemicaltransport model, VULCAN 2D, which lifts this restriction on wind patterns while sharing the same advantage as pseudo-2D models.While VULCAN 2D is constructed with horizontal and vertical dimensions in general, here we will focus on the zonal and vertical directions.
In Section 2, we describe the construction of the VUL-CAN 2D model.In Section 3, we validate the numerical scheme and modeling results through comparisons to an analytical solution, a pseudo-2D approach, and a 3D GCM, respectively.In Section 4, we delve into the applications to canonical hot Jupiters, HD 189733 b and HD 209458 b, including a super-solar C/O scenario for HD 209458 b.We use limiting cases to demonstrate the roles of horizontal and vertical transport and draw comparisons with previous works (Agúndez et al. 2014).We then compare VULCAN 2D to pseudo-2D models in Section 5, exploring when the uniform wind assumption in pseudo-2D models may no longer hold.Finally, we discuss the observational implications and the morning-evening limb asymmetry arising from horizontal transport in Section 6.

CONFIGURATION OF THE 2D CHEMICAL
TRANSPORT MODEL

The 2D Grid
A standard 1D photochemical kinetics model solves the continuity equation in the form of a set of coupled partial differential equations: where n is the number density (cm −3 ) of each species and t denotes the time.P and L are the chemical production and loss rates (cm −3 s −1 ) of the corresponding species at each vertical layer (Tsai et al. 2017).Equation ( 1) can be readily generalized to a 2D Cartesian space (x, z) to include horizontal transport: where ϕ z , ϕ x are the vertical (z) and horizontal (x) transport flux, respectively.In general, ϕ z encompasses vertical advection, eddy diffusion, and molecular diffusion (Tsai et al. 2021), whereas ϕ x describes horizontal transport on isobaric surfaces.Using isobaric flux is motivated by the fact that previous analysis of 3D GCMs is traditionally done in isobaric coordinates and winds on the isobaric levels typically dominate heat and chemical transport on tidallylocked planets (e.g.Koll & Komacek 2018;Showman et al. 2020;Hammond & Lewis 2021).We note that we use logpressure coordinates with the vertical coordinate defined as z = −H ln(p/P s ), where H is the pressure scale height that depends on the local temperature, mean molecular weight, and altitude-dependent gravity and P s is the reference pressure.While each vertical column uses the same pressure grid spaced uniformly in log-space, the corresponding z differs between columns due to the horizontal temperature gradient.Therefore, the horizontal derivative is evaluated at the same pressure level but not the same geometric height.The configuration of our 2D grid is illustrated in Figure 1.
The real challenge lies in numerically solving Equation (2).Most of the numerical methods for stiff equations require evaluating the Jacobian matrix (Brasseur & Jacob 2017).In a 1D system, the Jacobian matrix is neatly constructed through nested looping each species within each vertical layer (See Figure 14. in Tsai et al. (2017)).For example, the Rosenbrock method with a band matrix solver is employed in VUL-CAN.However, this structure breaks down when extending to 2D.To take advantage of the established numerical solver built for 1D systems, we apply an asynchronous integrator within each timestep (Small et al. 2013).Specifically, the right-hand side of Equation ( 2) is first discretized in x as where i, j denote the horizontal, vertical indices, respectively, and the +1/2 and -1/2 represent the interfaces enclosing that layer.Equation ( 3) is evaluated for each species and each vertical column.Since the last term (i.e. the horizontal transport flux) on the right-hand side of Equation (3) only contributes to the diagonal elements in the Jacobian matrix, Equation (3) has the same numerical structure as that in a 1D system (Equation ( 5) in Tsai et al. (2017)).This allows us to evaluate each column, including horizontal transport, using the existing 1D solver.Specifically, Equation (3) is computed for each x column within each timestep.Although errors can arise from the asynchronous update of horizontal transport flux associated with each column, we have tested integrating each column in different orders and found the errors to be negligible.
For the horizontal advection, we adopt a first-order upwind difference scheme, where the local concentration is affected by the upwind cell only (same as the vertical advection described in Tsai et al. (2021)).The advective flux for the k cell in the x direction is where v k−1/2 and v k+1/2 are the zonal wind velocity at the left and right interfaces of the grid cell k, respectively.In addition to advective transport, horizontal diffusion contributes to meridional transport on Jupiter and Saturn (Zhang et al. 2013;Hue et al. 2015;Hue et al. 2018) and is also implemented in VULCAN 2D.However, as horizontal winds can be directly obtained from GCM output, we focus only on advection for horizontal transport.For vertical transport, although the numerical functions include vertical advection, we restrict our model to vertical diffusion in this study and do not consider vertical advection to ensure numerical stability.We leave the full exploration of vertical advection to future work.
The advantage of our 2D spatial grid is it accommodates non-uniform two-dimensional flow patterns, such as the mean meridional or zonal circulation derived from the 3D GCM.In contrast, for pseudo-2D models that employ a Lagrangian rotating 1D column (Agúndez et al. 2014), the horizontal transport is restricted to a uniform jet by design.Both our 2D model and the pseudo-2D model share a limitation known as the Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1928), as already pointed out by Zhang et al. (2013).Given that we explicitly solve for each vertical column, the integration step size must be constrained by the time it takes for the horizontal flow to travel to neighbouring grid cells.This constraint becomes more stringent as the number of horizontal columns increases.In practice, to reduce the overall integration time to achieve convergence, we initially run our 2D model without horizontal transport (thereby removing the stepsize restriction) to attain the 1D steady state.We then run the full 2D VULCAN from this steady state (achieved without horizontal transport) with the stepsize limited by the CFL condition: dt < min( dx v x ) until the final steady state.Current VULCAN 2D evaluates each vertical column in sequence.Therefore, the simulation runtime scales with the number of columns linearly.However, there are opportunities to enhance code efficiency with parallel integration across the horizontal domain.Future development of VULCAN 2D will explore multiprocessing to speed up column-based computations.showing the volume mixing ratio at ξ = 1 with q = 0 (no zonal transport) and q = 1 (eastward zonal transport).The q = 1 case represents zonal transport on a timescale comparable to the chemical timescale.

Validation with analytical solutions
Zhang et al. ( 2013) derived an analytical solution for a 2D advective-diffusive system with parameterized chemical sources and sinks.We apply their analytical solution for a 2D zonal plane to validate our numerical model.To briefly recap, Equation (13) in Zhang et al. (2013) describes the steadystate abundance as a function of longitude and altitude, governed by vertical diffusion, horizontal advection, and chemical sources and sinks.We assume vertically uniform diffusion to simplify the expression, i.e., with α = γ = 0, the analytical solution for the volume mixing ratio under day-night sinusoidal chemical production (k = 1) is where χ is the volume mixing ratio, λ is longitude, ϕ is the phase shift, ξ is the vertical coordinate, P 0 and L 0 are the production and loss rates, respectively, and q relates the ratio of advection to chemical timescales, following the same notation as Zhang et al. (2013).
To compare with the analytical solution, we implemented a mock chemical network with only A and B as reactive species.The rate constant of A − −− → B is given by L 0 N 0 e −ξ /[A] and that of B − −− → A is given by P 0 N 0 (1 + coskλ)/[B], corresponding to the loss and production prescriptions in Section 5 of Zhang et al. (2013).We assume constant zonal winds of 71.5 m/s and R p = R J , together with P 0 N 0 = 2 × 10 −7 cm −3 s −1 and L 0 N 0 = 2 × 10 −7 cm −3 s −1 to yield q = 1.The vertical diffusion is set to 10 8 cm −2 s −1 but the horizontal distribution does not directly depend on this value.
Figure 2 compares our 2D numerical model to the analytical solution.The impact of the horizontal wind is twofold: it both transports and homogenizes the chemical gradient.While our numerical results show somewhat lower peak amplitudes than the analytical solution, which can likely be attributed to numerical diffusion, they effectively capture the overall distribution shape and correctly reproduce the "phase shift" due to eastward transport.

3.2.
Comparison with the Pseudo 2D approach (Agúndez et al. 2014) Agúndez et al. (2014) apply a rotating 1D model to represent an air column moving with a uniform zonal jet in a Lagrangian frame.Their chemical kinetics scheme was based on Venot et al. (2012).For a like-for-like comparison, we adopt the same rotating 1D-column as Agúndez et al. (2014) (referred to as Lagrangian-1D VULCAN) and compare it against our 2D photochemical model (2D VULCAN) using the identical N-C-H-O chemical scheme in VULCAN.We chose the canonical hot Jupiter HD 189733 b for this comparison, employing the same GCM output as that in Agúndez et al. (2014).In Lagrangian-1D VULCAN, the 1D column is set to travel around the equator over 2.435 days, following Agúndez et al. (2014), i.e. with a constant angular velocity of 1.493×10 −5 rad/s.To translate this rotation into the equivalent zonal-mean wind, we use a quasi-uniform zonal wind (faster at higher altitudes) in 2D VULCAN to match the same angular velocity.The zonal wind is 2430 m/s at 1 bar and scaled to each altitude according to the geometry.The equatorial plane is divided into 64 columns in longitude for both 2D VULCAN and Lagrangian-1D VULCAN.We used the same eddy diffusion coefficient (K zz ) for vertical mixing as Agúndez et al. (2014): K zz = 10 7 × (P/ 1 bar) −0.65 cm 2 s −1 .Other physical and chemical parameters were kept as much the same as possible between Lagrangian-1D and 2D VUL-CAN models.The comparisons between our Lagrangian-1D and 2D VULCAN for species displaying apparent longitudinal gradients are shown in Figure 3. 2D VULCAN and Lagrangian-1D VULCAN exhibit consistent results for species with either short (e.g.H) or long (e.g.CH 4 ) chemical timescales.This indicates that 2D VULCAN correctly captures the processes of vertical mixing and horizontal transport.We find minor deviations that appeared in trace species (volume mixing ratios ≲ 10 −6 ).For instance, NH 3 exhibits the most notable differences above 10 −4 bar.The discrepancy is likely due to the numerical integration: the first-order upwind in 2D VULCAN has stronger numerical diffusion than the Rosenberg method with second-order convergence in time used in Lagrangian-1D VULCAN.Nevertheless, the agreements are generally within a factor of two in the region below 0.1 mbar level.Our 2D VULCAN, with the equivalent uniform wind, can successfully reproduce the results obtained by the Lagrangian-1D approach.

Comparison with 3D transport (without photochemistry)
We have validated VULCAN 2D through an analytical solution and a pseudo-2D approach, to ensure the accuracy of our physical and numerical implementation within the 2D framework.In our next step, we will perform additional comparisons with a 3D GCM to see whether our 2D model effectively captures the primary transport process.
Mendonc ¸a et al. (2018b) previously investigated the global transport of chemically active tracers: H 2 O, CH 4 , CO, and CO 2 in the 3D simulation of WASP-43b.The chemical relaxation scheme (Tsai et al. 2018) (2018a,b).The eddy diffusion coefficients are often estimated from mixing length theory, using root-mean-square vertical velocity (w rms ) as the turbulent velocity and the scale height (H) as the mixing length.However, Parmentier et al. (2013) demonstrated that this approach tends to overestimate vertical mixing efficiency for hot Jupiters.Given this likely overestimation, we make two assumptions: K zz = w rms H and K zz = 0.01×w rms H to bracket the plausible range of K zz .The average temperature from the GCM is used to compute the scale height H = k B T mg (k B : Boltzmann constant; g: altitudedependent gravitational acceleration), but we fix the mean molecular weight m to 2.2387, following the same setting in Mendonc ¸a et al. (2018a,b).
Figure 4 compares the abundance distribution of CH 4 (for solar composition) and H 2 O (for C/O = 2) in the four quadrants of the equatorial region computed by VULCAN 2D and Mendonc ¸a et al. (2018a).Our 2D model with K zz = 0.01 × w rms H predicts vertical quench levels close to those in the 3D results.The preference for a weaker K zz might be attributed to the high gravity of WASP-43b.The slight increase of H 2 O between 0.1 and 0.01 in the 3D GCM is likely associated with meridional transport from higher latitude (Drummond et al. 2018) not accounted for in the 2D model.The overall good agreement between VULCAN 2D and the 3D results reinforces our 2D model's representation of the mixing process in the equatorial region of a hot Jupiter.

RESULTS: EQUATORIAL CHEMICAL DISTRIBUTIONS ON HD 189733 B AND HD 209458 B
We now first present an overview of the abundance distribution in the atmosphere of our fiducial simulations of HD 189733 b and HD 2094598 b.We briefly compare our results to previous work by Agúndez et al. (2014) on the same planets.To gain insight into the effects of vertical and horizontal transport, we examine the limiting cases where individual transport processes are isolated following the approach of Agúndez et al. (2014).We then explore the sensitivity to the eddy diffusion coefficients converted from the GCM wind and discuss the results guided by the limiting cases.Lastly, we explore the global chemical transport when C/O ≳ 1 with HD 209458 b, motivated by recent high-spectral-resolution observations.

Setup and overview
We adopt 3D GCM output from SPARC/MITgcm (Showman et al. 2009;Adcroft et al. 2004) for HD 189733 b and HD 209458 b.Chemical equilibrium is assumed when calculating radiative transfer with the correlated-k distribution method.For HD 189733 b, we use the same GCM output as Agúndez et al. (2014), where detailed GCM parameters are further listed in Table 1 of Steinrueck et al. (2019) (Table 1).For HD 209458 b, we do not include TiO and VO as shortwave opacity sources in our simulation, based on the lack of evidence of a thermal inversion from the reanalysis of emission spectra obtained by Spitzer and HST (Diamond-Lowe et al. 2014;Line et al. 2016).The parameters for our HD 209458 b GCM can be found in (Parmentier et al. 2013), with the change that TiO and VO opacities removed.The lower H and HCN, build up on the dayside while being transported horizontally by the zonal jet at the same time.We find that the chemical transport exhibits qualitative similarities between HD 189733b and HD209458 b, with the main difference being that CH 4 and NH 3 are more favored on the cooler HD 189733 b.
This behavior of chemical transport can be understood by comparing the horizontal and vertical dynamical timescales, as illustrated in Figure 11: Horizontal mixing dominates at altitudes below ∼0.1 mbar, while vertical mixing dominates at altitudes above ∼0.1 mbar level.Within the horizontal mixing region, species with longer timescales tend to display more uniform global abundances, whereas compositional gradients in the upper atmosphere are mainly controlled by vertical mixing and photochemistry.

Limiting cases with only horizontal transport and with only vertical mixing
Following the pedagogical exercise in Agúndez et al. (2014), this section presents the simulated vertical abundance profiles of HD 189733 b and HD 209458 b under limiting cases to isolate individual dynamical effects.Specifically, we examine scenarios where we exclude horizontal transport and vertical mixing, respectively.Figures 9 and 10 first show the distributions of CH 4 and H in these limiting cases, demonstrating the principles of horizontal transport and vertical mixing.As expected, vertical mixing tends to homogenize the vertical gradient of composition, while horizontal transport tends to homogenize the horizontal gradient.However, the influence of transport processes on a species depends on its chemical properties, which can be broadly categorized into two groups: For species like CH 4 that are replenished from the deep thermochemical regions, vertical mixing provides a zeroth order prediction on their abundances.For photochemical products like H that are predominantly produced on the dayside upper atmosphere, horizontal transport is crucial to account for their circulation to the nightside.
Detailed comparisons between the nominal 2D models, the vertical-mixing case, the horizontal-transport case, and chemical equilibrium at different longitudinal locations are presented in 12 and 13.Examining the quenched species such as CH 4 and NH 3 , we find that the vertical mixing and horizontal transport cases have the same quench levels.This is not surprising, as these quench levels correspond to the transition point associated with the same chemical timescales, regardless of the dynamical process.
However, the abundance distributions above the quench levels differ substantially between the vertical mixing and horizontal transport scenarios.This is because vertical mixing makes chemical abundances quenched from the hot deep layers, while horizontal transport results in quenching from the hot and irradiated dayside, as also noted in Agúndez et al. (2014).Owing to this nature of horizontal transport, the substellar abundance profiles closely resemble those from the vertical mixing case (the close match between the solid and dashed lines) for both planets.
Taking a closer look at Figures 12 and 13, species with long chemical timescales like CH 4 and NH 3 closely follow the vertical mixing distribution at hotter longitudes, the substellar point and evening limb, up to about 10 −4 bar.On the other hand, at the colder antistellar point and morning terminator, the equilibrium abundances differ substantially from the dayside, causing the abundance distribution predicted by the vertical-mixing model to diverge from the nominal 2D model (with both vertical and horizontal transport) at these cooler locations.These results shed light on the limitation of applying 1D models to interpret observations that probe the chemical properties of different regions.
One crucial feature of horizontal transport is transporting the photochemical products from the dayside to the nightside (Agúndez et al. 2014;Baeyens et al. 2022).For both planets, atomic H produced by photolysis can penetrate into the nightside, even reaching regions around the antistellar point where no UV photons are available.Transport of atomic H plays a key role in reacting with CH 4 and NH 3 to form HCN on the nightside.The abundances of CH 4 and NH 3 in Figures 12 and 13 fall between those in the vertical mixing and horizontal transport cases, highlighting the importance of considering both mixing processes.
To conclude our limiting-case analysis, we note that for hot Jupiters similar to HD 189733 b and HD209458 b, 1D model including vertical mixing serves as a fairly good approximation for the dayside.Given that the equatorial jet efficiently transports heat, 1D models generally better capture the hotter evening limb better compared to the cooler morning limb in the pressure regions most relevant for transmission spectroscopy.We will discuss the limb asymmetry as a result of horizontal transport more in Section 4.5.2 and 6.

Comparison with Agúndez et al. (2014)
There are several differences in modeling assumptions and setups between our models of HD 18933b and HD 209458 b and those presented in Agúndez et al. (2014): (i) Agúndez et al. (2014) assume uniform zonal winds, while we adopt longitude-and pressure-dependent wind profiles from the GCM.(ii) Shortwave opacity sources of TiO and VO, which are responsible for generating thermal inversion, are included in Agúndez et al. (2014) but are not included in our model for HD 209458 b (iii) Agúndez et al. (2014) use the chemical scheme of Venot et al. (2012), while our study adopts the chemical scheme of Tsai et al. (2021).Despite these differences, our nominal 2D VULCAN and those from the limiting cases are qualitatively consistent with the pseudo-2D outcomes presented in Agúndez et al. (2014).The major difference lies in the thermal inversion of HD 209458 b in Agúndez et al. (2014), making the dayside temperature at low pressures ∼ 1000 K higher than that in our model.Consequently, Agúndez et al. (2014) predict lower levels of CH 4 , NH 3 , and HCN on HD 209458 b.The hotter dayside also leads to a lower equilibrium abundance of CO 2 , creating the notable CO 2 day-night contrast seen in Agúndez et al. (2014).Overall, apart from the temperature inversion included in Agúndez et al. (2014), we show qualitative agreements in terms of quenching and transport behaviors.We will discuss the implications of assuming uniform zonal winds (i.e.pseudo-2D model) in Section 5 in more detail.

Sensitivity to K zz
While VULCAN has the capability to employ vertical advection instead of eddy diffusion (Tsai et al. 2021), it is numerically more stable to employ diffusion in the vertical direction in conjunction with molecular diffusion.One appar- ent caveat is that the parameterization of vertical mixing with eddy diffusion has been a long-standing uncertainty in atmo- spheric modeling (e.g., Smith 1998;Parmentier et al. 2013;Zhang & Showman 2018;Komacek et al. 2019).To test the sensitivity to uncertainties in estimating the eddy diffusion coefficient (K zz ), we vary the eddy diffusion coefficient profile by an order of magnitude in our model of HD 189733 b, both increasing and decreasing it.This explored range roughly spans the uncertain range considered in the previous literature (Smith 1998;Parmentier et al. 2013).
We examine the effects of varying vertical eddy diffusion coefficients on CH 4 and NH 3 , the major carbon and nitrogen species that exhibit transport-induced disequilibrium, in Figure 14.Stronger vertical mixing makes the transition to a vertical-mixing-dominated region occur at a higher pres-sures, as expected.This allows morning-evening asymmetries to emerge at somewhat higher pressures when horizontal transport becomes less efficient compared to vertical mixing.The pressure levels where morning-evening asymmetry becomes notable are generally between 1 and 0.05 mbar, as summarized in Table 1.Since the equilibrium abundances of CH 4 and NH 3 monotonically increase with pressure in our HD 189733 b model (also see the discussion in 4.3.2. in Tsai et al. (2021) and Fortney et al. (2020)), increased vertical mixing also results in slightly higher quenched abundances, while weaker mixing leads to lower abundances from horizontal quenching from the dayside.Other hydrocarbons and HCN also follow similar trends in response to mixing processes as CH 4 and NH 3 .Despite these effects, the global distribution of main species only mildly depends on the choice of K zz within the explored range.For HD 189733 b and HD 209458b, the bulk of the atmosphere below the millibar level still remains in the horizontal-transport-dominated regime, even when considering strong eddy diffusion.The transit observations using high-resolution spectroscopy of HD209458 b reported detection of H 2 O, CO, HCN, CH 4 , NH 3 , and C 2 H 2 molecules at statistically significant levels (Giacobbe et al. 2021).Notably, the presence of CH 4 , C 2 H 2 and HCN are indicative of super-solar carbonto-oxygen ratio (C/O), consistent with the previous highresolution analysis that also detected HCN (Hawker et al. 2018).In contrast, recent JWST/NIRCam transmission spectrum of this planet reported non-detections of CH 4 , HCN, and C 2 H 2 with upper limits provided (Xue et al. 2023), interpreted to have sub-solar C/O based on chemical equilibrium.While the interpretation by Giacobbe et al. ( 2021) is driven by cross-correlating with grids of models also assuming thermochemical equilibrium, CH 4 is influenced by quenching while C 2 H 2 and HCN can be efficiently produced by photochemistry out of equilibrium.In addition to the nominal HD 209458 b model with solar metallicity and C/O presented in this section, we run a 2D VULCAN model for HD 209458 b with super-solar C/O to explore the chemical transport under the alternative scenario of a carbon-rich atmosphere.We follow the best-fit C/O value of 1.05 as determined by Giacobbe et al. ( 2021) and keep the same solar metallicity for comparison with our nominal model, since metallicity is not well constrained in Giacobbe et al. (2021).

Enhanced hydrocarbon and horizontal transport induced limb asymmetry
The top and middle panels of Figure 15 depict the same mixing-ratio distributions as in Figure 8 but for C/O = 1.05.When C/O exceeds unity, H 2 O loses its dominance to CH 4 due to the lack of available oxygen after CO (Madhusudhan 2012;Moses et al. 2013;Heng et al. 2016).What is intriguing is the production of water resulting from the photolysis of CO in the dayside upper atmosphere above 0.1 mbar.The same process also operates in a solar C/O condition (Moses et al. 2011), but the increase of H 2 O is more pronounced here, owing to its lower equilibrium composition.Compared to solar C/O, C 2 H 2 abundances are significantly enhanced.Moreover, C 2 H 2 exhibits substantial zonal variations, spanning several orders of magnitude across the equator between 1 bar and 0.1 mbar level (bottom panel of Figure 15).
Although CH 4 is rather uniform below 10 −4 bar across the planet, several hydrocarbons show notable compositional gradients in the zonal direction, as seen in Figures 15 and  16.C 2 H 2 abundance on the evening limb is about 100-1000 times higher than that on the morning limb, whereas C 2 H 4 peaks between 0.1 and 10 −3 bar on the morning limb with strong vertical variations.The morning-evening limb asymmetry in C 2 H 2 is driven by the temperature difference and manifested due to its relatively short chemical timescale.C 2 H 2 is favored at the hotter evening limb, where the temperature is about 200-300 K higher than the morning limb, as evident in the profiles without zonal wind in Figure 16.In this region, the main destruction pathway for C 2 H 2 proceeds through a sequence of unsaturated hydrocarbons: The timescale of C 2 H 2 can be estimated from the ratelimiting steps, either the formation of CH 3 or C 2 H 3 in the above pathway.At 1 mbar pressure level, the lifetime of C 2 H 2 on the nightside ranges from 10 3 to 10 5 seconds.This timescale is comparable to that of horizontal transport (Figure 11), explaining the longitudinal compositional gradient C 2 H 2 displays.
Although C 2 H 4 also has a higher equilibrium abundance on the warming evening limb, photochemistry and horizontal transport lead to an accumulation of C 2 H 4 on the cooler morning limb instead.The combining effects result in a peak C 2 H 4 distribution on the morning limb greater than the abundance on the evening limb, in contrast with the lower morning C 2 H 4 abundance predicted in the absence of zonal transport (bottom panel of Figure 16).The photolysis of CH 4 on the dayside produced methyl radical (CH 3 ), the precursor to produce other hydrocarbons.CH 3 flow into the nightside where the lower temperature promotes the formation of C 2 H 6 , C 2 H 5 , and C 2 H 4 , since the combination of CH 3 into C 2 H 6 is exothermic and kinetically favored at cooler temperatures.With the aid of horizontal transport, the photochemically produced CH 3 is able to transport to the nightside to initiate the production of hydrocarbon species.This hydrocarbon production on the nightside is then carried to the morning limb by the zonal wind, leading to the peak shown in C 2 H 4 .Similar behavior is found in C 2 H 6 distribution as well, but at significantly lower abundances.For completeness, the equatorial abundance distribution as a function of the longitude and pressure of several key species can be found in Figure 22.

WHEN DOES THE ASSUMPTION OF A UNIFORM ZONAL JET IN PSEUDO-2D MODELS BREAK DOWN?
The equatorial jet is a robust feature for tidally-locked planets that receive steady day-night thermal forcing (Heng & Showman 2015;Showman et al. 2020).However, the equatorial jet transitions to a day-to-night flow when the radiative timescale or drag timescale becomes short (Showman et al. 2015;Komacek et al. 2019).In the case of cooler sub-Neptunes or nonsynchronously rotating planets, the zonal wind in the equatorial region can also develop a more complex structure with winds changing directions across varying pressure levels (Showman et al. 2015;Carone et al. 2020;Charnay et al. 2020;Innes & Pierrehumbert 2022).Here, we examine how horizontal transport changes on a hot Jupiter with strong frictional drag dominated by a day-to-night flow.We adopted the output of T eq = 1600 K with strong drag (τ = 10 4 s) from Tan & Komacek (2019) as our fiducial atmosphere with a day-to-night flow.The equatorial thermal and wind structures of our fiducial day-to-night circulation and those of HD 189733 b are compared in Figure 23.Our goal is to determine when a pseudo-2D approach with uniform flow remains valid and when it breaks down.

Comparisons within the superrotating regime
The equatorial jets induced by the stationary day-night heating in our HD 189733 b and HD 209458 b models resemble Gaussian-like vertical profiles, gradually diminishing towards zero in both the upper and deeper layers of the atmosphere, as depicted in Figure 5.We use HD 189733 b as an example of circulation characterized by equatorial superrotation for the comparison between 2D and pseudo-2D approaches.For the uniform zonal wind, the pseudo-2D model needs to adopt either the peak jet speed (Agúndez et al. 2014) or the averaged wind velocity over the jet region (Baeyens et al. 2021).For our pseudo-2D model of HD 189733b, we adopt a mean zonal wind of 1736 m/s at 1 bar from the GCM as the uniform zonal wind.
Figure 17 shows the comparison of CH 4 distributions, demonstrating differences caused by the simplified transport in the pseudo-2D approach verse the 2D model.In the deep region where thermochemical equilibrium holds, the actual wind pattern is not effectively relevant.The CH 4 distribution transitions from a horizontal homogenized regime to a vertical mixing-dominated regime at the same pressure level around 10 −4 bar in both models.In the upper atmosphere above 10 −4 bar, the pseudo-2D model has slightly slower horizontal transport (Figure 11) and predicts slightly more vertically mixed profiles compared to those in our nominal 2D model, where stronger winds around 1 mbar level at certain longitudes can alter the distribution.Similar trends are seen in other species as well.Despite these minor differences, we find the pseudo-2D approach to be a valid assumption when a broad superrotation jet is present.

comparisons within the day-to-night flow regime
Next, we explore the scenario of a hotter hot Jupiter with strong radiative (e.g., Komacek & Showman 2016) or magnetic drag (e.g., Perna et al. 2010), where the global circulation has changed from equatorial superrotation to day-tonight flow (Miller-Ricci Kempton & Rauscher 2012;Showman et al. 2015;Tan & Komacek 2019).Figure 18 compares the abundance distributions of CH 4 and HCN, the two main molecules that display compositional gradients, in the equatorial region of the fiducial strong-drag hot Jupiter atmosphere with T eq = 1600 K simulated by our nominal 2D model and pseudo-2D model.It is evident that the day-night flow leads to a symmetrical distribution, differing from the distribution governed by uniform zonal winds.This discrepancy is most pronounced around the morning limb, owing to the distinctive transport dynamics at play.In the nominal model, nightside-to-morning-limb advection occurs, contrasting with the morning-limb-to-nightside advection in the model with uniform zonal winds.This circulation pattern also enables more efficient transport of atomic H produced photochemically on the dayside to the morning limb.In the nominal 2D model, the CH 4 and HCN abundances remain low on both morning and evening limbs.Conversely, when assuming uniform eastward winds, CH 4 and HCN exhibit higher abundances around the morning limb due to nightside transport.Similar to CH 4 , species favored on the nightside over the dayside, such as NH 3 , follow similar trends showing higher morning limb abundances in the uniform winds model as well.
The disparity in abundance profiles across longitudes from the nominal model and pseudo-2D model is further illustrated in Figure 19.Except near the substellar point and at high pressures (P ≳ 1 bar), this comparison demonstrates that for species susceptible to photochemistry, the assumption of uniform zonal winds in the pseudo-2D framework can yield orders of magnitude differences in the observable region (approximately 0.1 to 10 −5 bar) when the dominant circulation is featured by a day-to-night flow.The pseudo-2D approach

SYNTHETIC SPECTRA (FOR THE EVENING AND MORNING LIMBS)
We present synthetic transmission spectra for the morning and evening limbs based on our 2D model results.The transmission spectra are computed using PLATON (Zhang et al. 2019(Zhang et al. , 2020)).including opacity sources of CH 4 , CO, CO 2 , C 2 H 2 , H 2 O, HCN, NH 3 , O 2 , NO, OH, C 2 H 4 , C 2 H 6 , H 2 CO, NO 2 , and collision-induced absorption (CIA) of H 2 -H 2 and H 2 -He.We will focus on the observational impact of horizontal transport and the differences between the evening and morning limbs.It is worth noting that the distribution of clouds and hazes can also significantly influence limb asymmetry (Kempton et al. 2017;Powell et al. 2019;Steinrueck et al. 2021;Savel et al. 2023).For the scope of this study, we will leave the effects of clouds to future work and will solely delve into the chemical transport within our cloud-free models.
Figure 20 illustrates the synthetic transmission spectra with and without horizontal transport for HD 189733 b and HD 209458 b (including solar and super-solar C/O).For HD 189733 b, 1D models (i.e.without horizontal transport) produce methane features on the cooler morning limb that are absent on the evening limb (Figure 12).However, this compositional gradient is readily homogenized once the zonal transport is included.For HD 209458 b (solar C/O), methane abundance remains too low in both 1D and 2D models (below ppm level; see Figure 13).Consequently, the influence of horizontal transport on the spectra of HD 209458 b with solar C/O is negligible.Instead, the dominant molecules that show up in the spectra are H 2 O, CO, and CO 2 , all of which exhibit relatively uniform equilibrium abundances throughout the planet.As a result, these gases do not contribute to limb asymmetry, whether horizontal transport is in play or not.In the case of HD 209458 b with super-solar C/O, CH 4 took over H 2 O to make the strongest spectral features at 2.3, 3.3-3.9,5.5-6.6, and 7-8.5 µm.As CH 4 becomes the predominant carbon-bearing molecule, it also reaches a uniform distribution across the entire planet.The major limb asymmetry is the presence of C 2 H 4 absorption on the morning limb but not on the evening limb, due to the horizontal transport of C 2 H 4 from the nightside to the morning limb.

DISCUSSIONS AND CONCLUSIONS
In this paper, we present the 2D version of the photochemical model VULCAN.We first validate VULCAN 2D with analytical solutions.VULCAN 2D successfully reproduces the special case equivalent to the pseudo-2D approach with uniform winds and also demonstrates consistent results with a 3D GCM (Mendonc ¸a et al. 2018b).We use limiting cases to demonstrate the distinct effects of the vertical and horizontal mixing processes.For canonical hot Jupiters, such as HD 189733 b and HD 209458 b, we find most of the atmo-sphere below 1-0.1 mbar is within the horizontal transportdominated region, where zonal advection prevails over vertical mixing.In the upper atmosphere above this region, photochemistry and vertical mixing control the composition.We explore the sensitivity to the parametrization of vertical mixing and find a mild dependence in the abundance distribution of our hot Jupiter models.We note that stronger vertical mixing can, in principle, promote morning-evening asymmetry.
For HD 189733 b, the morning-evening limb asymmetry in CH 4 predicted by 1D models is readily homogenized when horizontal transport is included.For HD 209458 b with solar C/O, the transmission spectra exhibit no limb asymmetry attributed to the composition due to the paucity of CH 4 .However, with super-solar C/O, horizontal transport results in notable limb asymmetries in hydrocarbons (C 2 H 4 in this case).For atmospheres with circulation dominated by an equatorial jet, we show that the pseudo-2D (rotation 1D column) approach can reasonably capture the transport, but the assumption of uniform flow breaks down for day-to-night circulations under stronger drag and pseudo-2D models can yield observable differences (0.1 to 10 −5 bar) orders of magnitude apart, particularly for species susceptible to photochemistry or with inherent compositional gradient in equilibrium.
The 2D modeling framework highlights the need to consider both horizontal and vertical transport when interpreting the compositions from transmission observations probing the limbs.The 2D framework developed here bridges the gap between traditional 1D photochemical kinetics models and 3D general circulation models that typically exclude chem-ical kinetics.Future directions include incorporating sulfur chemistry and applying the model to the meridional plane or tidally locked coordinate (Koll & Abbot 2015) to explore the role of overturning circulation.
In this work, we did not explore the effects of vertical advection.The upward or downward transport can lead to significantly different distributions, as already demonstrated in the ammonia distribution on Jupiter by a 1D model (Tsai et al. 2021).Hammond & Lewis (2021) also showed how slow overturning winds could transport much more heat than fast zonal jets on tidally locked planets with weak temperature gradients.A future model development should represent both of these processes.
By elucidating how the atmosphere composition is regulated by global circulation, this 2D modeling approach will pave the way for self-consistent and more comprehensive models and provide a useful tool to enhance the capacity of 3D GCMs for interpreting observations.

Figure 1 .Figure 2 .
Figure1.The grid configuration of VULCAN 2D (left) and the stencil of the upwind scheme for horizontal (i) advection and vertical ( j) diffusion (right).The arrows illustrate the advection and diffusion transport between adjacent grids.

Figure 3 .
Figure3.Mixing ratio profiles at various longitudes of HD 189733 b using Lagrangian-1D VULCAN (dashed) compared to those using 2D VULCAN with equivalent zonal winds (solid; see text for details).The substellar point is at 0 • longitude.Lagrangian-1D VULCAN has the same configuration as the pseudo-2D model (e.g.,Agúndez et al. 2014;Baeyens et al. 2021).Equilibrium abundances are indicated by the grey lines.
implemented in Mendonc ¸a et al. (2018b) simplifies the thermochemical kinetics with a linear response according to the chemical timescale.We compare with the work by Mendonc ¸a et al. (2018b) because their chemical timescales are derived from the same chemical network implemented in VULCAN.Here, photochemistry in 2D VULCAN is switched off to have an equivalent compari-

Figure 5 .Figure 6 .Figure 7 .
Figure 5.The equatorial temperature (left) and zonal winds (right) from our GCMs of HD 209458 b and HD 189733b, averaged over ± 30 • latitudes.The substellar point is located at 0 • longitude.The uniform zonal winds (constant angular velocities with 2430 m/s for HD 189733 b and 3850 m/s for HD 209458 b at the 1 bar pressure level, respectively) assumed in Agúndez et al. (2014) are illustrated in black dotted lines in the right panels for comparison.and upper pressure boundaries of the GCM simulations are 200 bar and 2× 10 −6 bar, respectively.To extend the model domain for VULCAN 2D, the temperatures and winds at the upper GCM boundary are held constant up to 10 −8 bar.The equatorial average temperatures and winds for our HD 189733b and HD 209458 b simulations are shown in Figure5.The zonal wind speeds evidently vary across longitudes at pressures above 0.1 bar.Both planets share similar thermal and dynamical structures in the equatorial region, though HD 209458 b exhibits higher temperatures and overall faster zonal wind speeds.We adopt the same global eddy diffusion coefficients asAgúndez et al. (2014) for vertical mixing.These K zz profiles are compared with those estimated by the mixing length theory in Figure21.We will further explore the sensitivity to the adopted eddy diffusion coefficients

Figure 8 .
Figure 8.The vertical profiles in the equatorial regions of HD 189733b and HD 209458 b simulated by 2D VULCAN.The shades of colors from light to dark correspond to longitudes from the substellar point (0 • ) through nightside and back around to the substellar point (360 • ).

Figure 9 .
Figure 9.The equatorial distribution of CH 4 as a function of longitude and pressure on HD 189733 b, from our nominal model, vertical-mixing model, horizontal-transport model, and chemical equilibrium.Note the colorbars have different scales.

Figure 11 .
Figure 10.Same as Figure 9 but for H.

Figure 13 .
Figure13.The abundance profiles of selected species on HD 209458 b at four longitudes: substellar point, evening limb, antistellar point, and morning limb.The volume mixing ratios computed by the nominal 2D model (solid) are compared to the vertical-mixing case (excluding horizontal transport; dashed), the horizontal-transport case (excluding vertical diffusion; dash-dotted), and thermochemical equilibrium (dotted).

Figure 14 .Figure 15 .
Figure 14.CH 4 and NH 3 vertical profiles in the equatorial regions of HD 189733 b, simulated with 10 times stronger (orange) and weaker (green) eddy diffusion.The shades of colors from light to dark correspond to longitudes from the substellar point (0 • ) through nightside and back around to the substellar point (360 • ).4.5.Super-solar C/O of HD 209458 b 4.5.1.background

Figure 16 .
Figure 16.The vertical profiles of CH 4 , C 2 H 2 , and C 2 H 4 on the morning and evening limbs in our 2D model of HD 209458 b with C/O = 1.05.The distributions without zonal wind are shown in dashed lines, while chemical equilibrium abundances are indicated by dotted lines.

Figure 17 .
Figure17.The mixing ratio profiles of CH 4 across different longitudes computed from our nominal HD 189733 b model with pressure-and longitude-dependent zonal winds (solid) and from the pseudo-2D model with uniform zonal winds (dashed).We adopt a mean zonal wind velocity of 1736 m/s at 1 bar in the nominal model for the uniform wind speed in the pseudo-2D model, following the same approach as inAgúndez et al. (2014).

Figure 20 .
Figure 20.The synthetic transmission spectra of the evening and morning limbs of HD 189733 b and HD209458 b (for both solar and super-solar C/O).The atmospheric structure is derived from MITgcm simulation (Parmentier et al. 2013) with the chemical composition computed by 2D VULCAN.The models excluding zonal transport are also shown for comparison.

Table 1 .
The pressure levels (mbar) where the morning and evening abundances begin to deviate by more than 50%.
Figure12.The abundance profiles of selected species on HD 189733 b at four longitudes: substellar point, evening limb, antistellar point, and morning limb.The volume mixing ratios computed by the nominal 2D model (solid) are compared to the vertical-mixing case (excluding horizontal transport; dashed), the horizontal-transport case (excluding vertical diffusion; dash-dotted), and thermochemical equilibrium (dotted).