  Increasing main sequence stellar luminosity with stellar mass leads to the
eventual dominance of radiation pressure in stellar envelope hydrostatic
balance. As the luminosity approaches the Eddington limit, additional
instabilities (beyond conventional convection) can occur. These instabilities
readily manifest in the outer envelopes of OB stars, where the opacity increase
associated with iron yields density and gas pressure inversions in 1D models.
Additionally, recent photometric surveys (e.g. TESS) have detected excess
broadband low frequency variability in power spectra of OB star lightcurves,
called stochastic low frequency variability (SLFV). This motivates our novel 3D
Athena++ radiation hydrodynamical (RHD) simulations of two 35$\,$M$_\odot$ star
envelopes (the outer $\approx$15$\%$ of the stellar radial extent), one on the
zero-age main sequence and the other in the middle of the main sequence. Both
models exhibit turbulent motion far above and below the conventional iron
opacity peak convection zone (FeCZ), obliterating any ``quiet" part of the
near-surface region and leading to velocities at the photosphere of
10-100$\,$km$\,$s$^{-1}$, directly agreeing with spectroscopic data. Surface
turbulence also produces SLFV in model lightcurves with amplitudes and
power-law slopes that are strikingly similar to those of observed stars. The
characteristic frequencies associated with SLFV in our models are comparable to
the thermal time in the FeCZ ($\approx$3-7$\,$days$^{-1}$). These simulations,
which have no free parameters, are directly validated by observations and,
though more models are needed, we remain optimistic that 3D RHD models of main
sequence O star envelopes exhibit SLFV originating from the FeCZ.
