# Regulating star formation in a magnetized disk galaxy

Hector Robinson,1 James Wadsley,1
1Department of Physics and Astronomy, McMaster University, 1280 Main Street
West, Hamilton, Ontario, L8S 4M1 Canada

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use high-resolution MHD simulations of isolated disk galaxies to
investigate the co-evolution of magnetic fields with a self-regulated, star-
forming interstellar medium (ISM). The simulations are conducted using the
Ramses AMR code on the standard Agora initial condition, with gas cooling,
star formation and feedback. We run galaxies with a variety of initial
magnetic field strengths. The fields grow rapidly and achieve approximate
saturation within 500 Myr, but at different levels. The galaxies reach a
quasi-steady state, with slowly declining star formation due to both gas
consumption and increases in the field strength at intermediate ISM densities.
We connect this behaviour to differences in the gas properties and overall
structure of the galaxies. In particular, strong fields limit feedback
bubbles. Different cases support the ISM using varying combinations of
magnetic pressure, turbulence and thermal energy. Magnetic support is closely
linked to stellar feedback in the case of initially weak fields but not for
initially strong fields. The spatial distribution of these supports is also
different in each case, and this is reflected in the stability of the gas
disk. We relate this back to the overall distribution of star formation in
each case. We conclude that a weak initial field can grow to produce a
realistic model of a local disk galaxy, but starting with typical field
strengths will not.

###### keywords:

Methods: numerical – MHD – ISM: magnetic fields – Galaxies: star formation

††pubyear: 2023††pagerange: Regulating star formation in a magnetized disk
galaxy–Regulating star formation in a magnetized disk galaxy

##  1 Introduction

Magnetic fields have been detected at all scales in astrophysics (Han, 2017)
and are predicted to play important roles in galaxy evolution. A key question
is what role magnetic fields play in regulating galactic star formation. The
most straightforward consideration is that magnetic fields provide an
additional pressure which can support gas against gravitational collapse on
the scale of a galactic disk, but they can also affect galaxy-scale dynamics,
the properties of turbulence, the effectiveness of stellar feedback, and the
formation of molecular clouds. Importantly, turbulence acts to amplify and
reshape the magnetic field, so in practice we must study the joint evolution
of the magnetic field and the ISM together.

Detecting magnetic fields in galaxies is challenging. A diverse set of methods
allows us to infer differing aspects of the field. These include synchrotron
emission, the Zeeman effect, Faraday rotation, polarised thermal emission of
magnetically aligned dust grains (Pattle et al., 2023), and polarised emission
of starlight due to extinction by aligned dust grains.

Of primary interest is the field strength. Synchrotron observations can be
used to estimate total magnetic field strengths by assuming energy
equipartition between magnetic fields and cosmic ray particles. Fields
strengths measured in spiral galaxies with this method tend to be around
$\sim$10 $\mu$G and decrease slowly with galactocentric radius (Fletcher et
al., 2011; Basu & Roy, 2013; Beck & Wielebinski, 2013). The dependence on the
poorly constrained cosmic ray distribution means both the field strength and
gradient are highly uncertain. Due to the diffusive nature and resulting large
scale height of cosmic rays, it is expected to probe a thick volume around the
galactic disk (Zweibel, 2017).

In the Milky-Way and a few nearby galaxies, line-of-sight magnetic field
strengths have also been measured via the Zeeman Effect (Crutcher et al.,
2010), which causes emission lines to split when molecules are in the presence
of line of-sight magnetic fields. Zeeman observations in the Milky-Way have
found fields strengths of $\sim$ 10 $\mu$G at number densities of 10-100 cm-3,
In gas at number densities $\gtrsim$ 1000 cm-3 the upper envelope to the field
strength scales with number density as $B\propto n^{0.5-0.7}$ (Crutcher et
al., 2010). There is considerable spread in these measurements for 10-1000
cm-3. It is commonly assumed that typical field strengths remain flat at lower
densities.

The remaining techniques mostly indicate the field morphology. Many
observations find large-scale spiral patterns in galactic magnetic fields,
even in galaxies that do not have optical spiral structures (Chyży & Buta,
2008; Beck et al., 2019; Lopez-Rodriguez et al., 2023). Down to $\sim 100$ pc
scales, the fields tend to be aligned parallel with the structure, such as low
density filaments (Goldsmith et al., 2008; Sugitani et al., 2011). These
results have been confirmed more recently by synchrotron maps of molecular
clouds in the Milky way (Planck Collaboration et al., 2016).

There are many theoretical predictions for how magnetic fields should evolve
and affect their host galaxies. Tiny cosmic seed fields are amplified
exponentially to detectable levels within a few Gyr (Geach et al., 2023). On
small galactic scales, the turbulent dynamo is expected be ubiquitous. It can
exponentially amplify field strengths over timescales of $\lesssim 10$ Myr,
saturating at level that is a fraction of the turbulent energy (Federrath et
al., 2011; Rieder & Teyssier, 2016). The $\alpha-\Omega$ dynamo is expected to
order the small scale turbulent fields into disk scale regular fields over Gyr
timescales (Brandenburg & Subramanian, 2005). At intermediate scales there is
also the gravitational-instability dynamo associated with spiral structures
(Riols & Latter, 2019). Amplification rates in simulations are still highly
dependent on numerical resolution and feedback methods (Rieder & Teyssier,
2016), which can make comparisons between different simulations difficult,
however the level they saturate at appears to be independent of those effects.
Thus saturated fields present an appealing target for study that is less
dependent on numerical method differences.

After saturation, the energy density of magnetic fields relative to other
sources can vary depending on the phase of gas. The diffuse medium can have
significant magnetic support, meaning that the magnetic pressure is comparable
to thermal pressure (plasma $\beta$ $~{}\sim$ 1). At higher densities, where
gas is colder, the primary support is turbulent and the clouds are typically
magnetically supercritical (E${}_{\textrm{mag}}$ < E${}_{\textrm{grav}}$
$\sim$ E${}_{\textrm{turb}}$) (see review by Beck & Wielebinski 2013).

Magnetic fields play are expected to play key role in how gas transitions
between phases within the ISM (Krumholz & Federrath, 2019). This manifests as
a difference in the distribution of gas densities in the ISM, which are
created by turbulent compression and expansion. Turbulence predicts a
lognormal probability density function (PDF) that is also seen in observations
(Kainulainen et al., 2009). At high densities when the gas becomes
gravitationally unstable, it diverges from the lognormal (Burkhart, 2018).
Magnetic pressure narrows the width of the PDF by resisting turbulence’s
ability to compress gas. Simulations on cloud and kpc scales have shown that
this effect can reduce star formation rates by a factor of 2-3 (Federrath &
Klessen, 2012; Padoan et al., 2012; Girichidis et al., 2018; Krumholz &
Federrath, 2019; Kim et al., 2023). Turbulence also plays a role in the
magnetic field strength vs. gas density scaling relations. The 0.5-0.7 power
law was originally thought to come from gravitational contraction, but cloud-
scale simulations by Cao & Li (2023) retrieve the same scaling without self
gravity and the authors argue it comes from turbulent compression instead. A
key takeaway from this discussion is that the role of magnetic fields on
smaller scales is complex and under intense study. Working on slightly larger,
galactic scales, allows for a simpler treatment of star formation.

MHD simulations are a powerful tool for studying magnetic fields on many
scales. On cosmological scales, magnetic fields are generally not dynamically
dominant but some simulations have begun to include them (Grand et al., 2017;
Hopkins et al., 2018; Steinwandel et al., 2022; Martin-Alvarez et al., 2018).
Cosmological seed fields are the origin of all galactic fields; thus
cosmological simulations can be used to test models of magnetic field
amplification and produce toroidally dominated fields similar to those seen in
observations (Rieder & Teyssier, 2017b). One difficulty with cosmological
simulations is that magnetic fields amplify strongly during the poorly
resolved and highly chaotic infall and merging phases. These processes apply
ongoing, significant perturbations to the state of the galaxy. There is also
the computational expense of simulating a large cosmological environment for
the full history of the universe.

An alternative approach is high resolution simulations of individual galaxies,
including the effects of magnetic fields. Körtgen et al. (2019) simulated an
isolated disk galaxy and showed that magnetic fields can speed up disk
fragmentation and drive outflows hundreds of parsecs above the disk even
without stellar feedback. Galaxies that have both star formation and feedback
tend to have lower star formation rates when MHD is also included and can
magnetize their CGM with magnetic outflows (Steinwandel et al., 2019, 2020;
Pakmor et al., 2020; Wissing & Shen, 2023), suggesting that MHD can help
regulate star formation on disk scales. Recent simulations also suggest that
the magnetic field strength continuously declines with density, effectively as
a power law, extending to densities below the point where observationally
inspired models suggest it should become constant (e.g. (Ponnada et al.,
2022)).

Prior work has tended to focus on amplification rates and the final magnetic
configuration (e.g. Rieder & Teyssier, 2017a; Su et al., 2018). However, it is
also important consider how magnetic fields affect the state of the gas in the
multiphase interstellar medium and whether the star formation is regulated in
the same way as we introduce progressively stronger fields.

In this paper we present a controlled study, simulating an isolated galaxy
with a well-known initial condition, with and without magnetic fields, with
several different initial field strengths. Thus we can focus on the
development and co-evolution of magnetic fields due to self-regulated star
formation and feedback in a galactic disk. By running such cases for several
dynamical times, we expect to produce a steady-state, self regulated star-
forming ISM to study. In addition, by using a standard setup and simple, well-
tested star formation and feedback models, we aim to make the interpretation
more straightforward.

The remainder of the paper is organized as follows: In section 2, we describe
the simulation method and magnetized galaxy setup. In section 3.1, we analyse
the evolution of the magnetic field in each galaxy, including the approach to
a saturated state and comparing the strengths to observations. We examine the
resulting visual appearance of each galaxy in section 3.2. In section 3.3, we
examine the overall star formation and its radial distribution in each case.
We then explore how this is reflected in their ISM. Section 3.4 examines the
gas properties, seeking to determine the underlying drivers of the differences
in star formation rates and their connection to the magnetic fields and other
forms of support for the gas. In section 3.5, we study the combined effect of
the different support mechanisms on the gravitational stability of the
galactic disks. In section 4 we discuss our results and future work. Finally
we summarize our conclusions in section 5.

##  2 Simulation Method

We conduct magnetohydrodynamic (MHD) simulations of isolated galaxies using
the adaptive mesh refinement (AMR) code Ramses (Teyssier, 2002) to solve the
ideal MHD equations using an HLLD approximate Riemann solver (Miyoshi &
Kusano, 2005). The solenoidal constraint ($\nabla\cdot B=0$) is enforced with
the constrained transport method (Evans & Hawley, 1988). The dynamics of stars
and dark matter are solved using the particle-mesh technique (Hockney &
Eastwood, 1981). Gas cooling and heating is included via the Grackle chemistry
and cooling library (Smith et al., 2017). Grackle uses metal cooling rates
tabulated from output from the photo-ionization code Cloudy (Ferland et al.,
2017). We also include a photoelectric heating rate of $\zeta=4\times
10^{-26}$ erg cm-3 s-1, which allows for a two-phase ISM similar to that
proposed by Wolfire et al. (2003).

###  2.1 Initial Conditions and Refinement

The galaxies we simulate are all based on the medium-resolution isolated disk
galaxy from the Agora Project (Kim et al., 2016), but with an initial magnetic
field added. It has an active dark matter halo that follows an NFW profile, a
stellar bulge that follows a Hernquist profile, and a disk that is 80% stars
and 20% gas by mass with a density profile given by

$\rho_{\textrm{gas}}(r,z)=\rho_{0}e^{(-r/r_{d})}e^{(-|z|/z_{d})}$ (1)

With $\rho_{0}=M_{\textrm{gas}}/(4\pi r_{d}^{2}z_{d})$, where $r_{d}=3.432$
kpc and $z_{d}=0.1r_{d}$. The stellar and dark matter components are modelled
with collisionless particles, and the gas is initiated on a Ramses AMR Grid.
Kim et al. (2016) contains fulls details about the disk setup.

This initial condition has been described as similar to ’Milky-Way like’
spiral galaxy with a redshift of z$\sim$1\. However, the dark matter halo and
rotation curve are quite similar to a $z\sim 0$ large disk galaxy, such as NGC
5055 (the Sunflower galaxy). The initial surface densities of gas and stars
are also quite similar to NGC 5055 (as noted by Benincasa et al. 2020). In
this sense it is actually a reasonable proxy for a nearby spiral galaxy.

The gas disk is initialized inside of a domain 600 kpc on a side that has a
base grid of 643 cells, that is allowed to refine an additional 10 levels
which gives a spatial resolution of 9.15 pc at the highest level. A cell will
refine if it contains more than 10,000 M⊙ of gas or if it contains more than 8
collisionless particles. The disk starts off with $M_{\textrm{gas}}=$
8.59$\times$109 M⊙, which will decrease throughout the simulation as gas is
converted into stars. The gas is initialized to a temperature of 10,000 K and
solar metallicity.

Inside the gas disk, we initialize a magnetic field with a morphology that is
purely toroidal, containing no vertical or radial components. The field
strength scales with gas density as

$B=B_{0}\left(\frac{\rho}{\rho_{0}}\right)^{2/3}$ (2)

B0 is the value of the magnetic field strength on the midplane at the center
of the galaxy. The magnetic field strength decreases further out in the disk,
decreasing by a factor of 10 by 12 kpc. Due to the grid initialization
process, some initial B values are changed by $\pm$ 10%. Figure 1 includes
radial profiles of the initial magnetic field strength (straight red lines).

We simulate four galaxies which are identical except for the initial magnetic
field, with values for each case summarized in Table 1. The first case has
zero initial magnetic field, equivalent to being simulated with regular
hydrodynamics. The remaining three galaxies are referred to as MHD Weak, MHD
Medium, and MHD Strong. In each of them, the field strength initially scales
with gas density according to equation 2, but the constant $B_{0}$ is modified
so that the MHD Medium has fields 10 times stronger than the MHD Weak, and MHD
Strong has fields 10 times stronger than MHD Medium. An increase of 10 times
in magnetic field strength corresponds to an increase of 100 times in magnetic
energy.

Name | B${}_{\textrm{ISM}}$ | $\beta_{\textrm{ISM}}$ | B0
---|---|---|---
Hydro | 0 | $\infty$ | 0
MHD Weak | 0.1 $\mu$G | 800 | 0.85 $\mu$G
MHD Medium | 1 $\mu$G | 8 | 8.5 $\mu$G
MHD Strong | 10 $\mu$G | 0.08 | 85 $\mu$G
Table 1: Summary of initial magnetic properties of each simulation.
B${}_{\textrm{ISM}}$ and $\beta_{\textrm{ISM}}$ are magnetic field and plasma
$\beta$ values, respectively, in the typical ISM (gas density $n\sim 0.25$
cm-3). B0 is the corresponding magnetic field strength in equation 2 (at the
geometric center of the galaxy). Otherwise initial magnetic field strengths
scale as B$\ \propto\rho^{2/3}$. As a result of this scaling, the plasma
$\beta$ (P${}_{\textrm{thermal}}$/P${}_{\textrm{mag}}$) increases with radius.

###  2.2 Star Formation and Feedback

Stars particles are formed stochastically via a Schmidt-law of the form

$\frac{d\rho_{*}}{dt}=\frac{\epsilon_{\textrm{ff}}}{\rho\,t_{\textrm{ff}}}\quad\text{if}\quad\rho>\rho_{\textrm{crit}}$
(3)

where $\rho_{*}$ is the stellar density, $\epsilon_{\textrm{ff}}$ is the star
formation efficiency per free-fall time which we set to $\epsilon_{ff}=0.1$,
and $\rho_{\textrm{crit}}$ is a threshold density which corresponds to a
number density of 100 cm-3.

Stellar feedback is entirely supernovae (SN), injected as thermal energy 5 Myr
after stars first form with $10^{51}$ erg per 91 M⊙ of young stars. This
energy is treated via the delayed cooling model of Agertz et al. (2011), which
allows unresolved superbubbles to grow correctly by initially treating hot
supernova ejecta as unresolved, non-cooling bubbles whose energy is converted
to regular thermal energy with a e-folding time of 5 Myr. At the chosen
resolution, other forms of feedback are largely unresolved, as are the dense
structures in which they would chiefly operate.

We note that these choices differ from the original Agora simulations (Kim et
al., 2016). In particular, Agora used a low efficiency ($\epsilon_{ff}=0.01$)
in dense gas. This pushes the characteristic time for star formation to $\sim
1$ Gyr, effectively forcing a typical overall galactic star formation rate
(SFR), that was the same with or without feedback. With $\epsilon_{ff}=0.1$,
stellar feedback is necessary to regulate star formation to the level expected
for typical disk galaxies (Robinson, 2021). Feedback is expected to couple in
such a way to reproduce large scale ISM properties including the scale height
(Ostriker et al., 2010; Benincasa et al., 2016). Thus we have adopted a robust
and easy to reproduce approach to star formation and feedback for this work.

Each galaxy was evolved for 1 Gyr, corresponding to a handful dynamical times
at the outer radius. The inner regions were expected to evolve rapidly and
then settle into a quasi-steady state of ongoing star formation. Thus the goal
was to produce an interval of several 100 Myr to study this relatively quiet,
self-regulated state in each case.

##  3 Simulation results

At the start of the simulation, all four galaxies begin their evolution by
settling from the slightly unstable intial state by compressing vertically.
While this happens cold-phase gas condenses and fragments along spiral arms
into individual clouds. In the Hydro and MHD Weak cases, this collapse is
violent and results in a starburst, mostly localized to the galactic center (r
< 2 kpc). In the cases with stronger fields, the fields resist the
compression, preventing a starburst from occurring, and the onset of star
formation is delayed. As a result, the Hydro and MHD Weak galaxies have a
smaller fraction of gas remaining in the galactic center for the remainder of
the simulation, for this reason we exclude the center 2 kpc from some of our
analysis to ensure a fair comparison between the galaxies.

###  3.1 Magnetic Field Evolution

Figure 1: Magnetic field strength vs. galactocentric radius in each MHD
galaxy. Top row shows a mass-weighted average, and bottom row shows a volume-
weighted average. Color shows the the time of the snapshot from 0 (red), to 1
Gyr (pink). Black lines are data from nearby galaxies from Basu & Roy (2013).

The primary focus of this work is how magnetic fields influence the structure
of the ISM and the regulation of star formation at later stages. Before
examining that we give a brief overview of the magnetic field evolution
leading to the final state.

Figure 1 shows the average magnetic field strength vs. radius in each galaxy
over time, and compares them to field strengths from synchrotron observations
by Basu & Roy (2013). Firstly, we note that this observational sample includes
NGC 5055, for which the simulation set-up is a reasonable proxy. Secondly, the
variation among the observational sample is also quite small, with estimated
field strengths of $\sim$ 10-20 $\mu$G in the radii of interest, 2-10 kpc.

The figure includes two methods of calculating the average field strength. The
first is a mass-weighted average (top panels) and the second is a volume
weighted average (lower panels). The method of averaging does not particularly
affect the field estimate in the MHD Strong galaxy because the gas scale
height is quite large (see 3.4 for a quantitative comparison). For weaker
fields, the galaxies are thinner and the differences more pronounced.

The volume-weighted average in the MHD Weak galaxy is the only one that
saturates with a flat profile, as seen in the sychrotron estimates, albeit
with a lower field strength. Synchrotron emission arises from cosmic rays
which are typically assumed to have a large scale height (Zweibel, 2017),
making the volume-weighted average the closest match. Precise cosmic ray
energy distributions are hard to pin down. Regardless of the assumed cosmic
ray to magnetic energy ratio, the observed flat radial distribution should be
approximately preserved. This makes the MHD Weak case the most compelling. In
particular, we infer that initial fields must be well below equipartition in
order to naturally evolve to a realistic saturated state.

Similar amounts of magnetic field amplification are visible in both the MHD
Weak and MHD Medium galaxies regardless of the choice of weighting. The MHD
Strong galaxy experiences a net loss of field strength over time, confirming
that it was oversaturated from the beginning. Both the MHD Weak and Medium
galaxies appear close saturation after about 500 Myr, but they do not saturate
at the same value. Field strengths may decrease over time due to magnetic flux
leaving the disk. In the inner regions of the MHD Weak galaxy, field strengths
peak at around t=500 Myr and then decrease slightly due to flux leaving in
vertical outflows. In the MHD Strong galaxy, net flux loss is expected due to
magnetic braking.

In all cases, the amplification rate is lower in the outer disk, where there
is less star formation (see section 3.3). This behaviour is most consistent
with a turbulent style dynamo, both in terms of the rate of growth and the
strong association with feedback from star formation.

Our field strengths can also be compared to those inferred from Zeeman
measurements, by placing them on a field strength vs. density (B vs. n) plot,
shown in figure 2. Here we plot the median field strength in each density bin,
but variations of an order of magnitude are common at a given density. In this
figure, field amplification due to dynamo action is seen as a vertical
translation. When field strengths increase due to gas compression (or
expansion), the points should also move to the right as the density increases
(or decreases, respectively). This largely explains the rapid expansion of the
plot upward in density and field strength after the initial state. In our
galaxies, it is clear that most amplification is happening in the diffuse
medium with number densities of $0.1-10$ cm-3, while the dense gas strengths
remain constant or even slightly decrease over time. The result is a flat
power law, close to 0.5 in the MHD Weak case. There is a trend to shallower
power laws as we progress from the MHD Weak to the MHD Strong case. Thus, only
the MHD Weak case displays the steep power law at higher densities inferred by
Crutcher et al. (2010).

In the MHD Strong galaxy, any gas that got more dense than the initial
condition has typically experienced field increases of only a factor of two at
most. This can be explained if the gas flows are directed mostly along field
lines, which is expected in regions which are magnetically dominated. This
also explains why the mass-weighted average gives similar field strengths in
figure 1, because there is hardly any amplification in the high-density gas.

At low densities, the simulated field strengths continue to decrease with
density, showing no hint of the constant field extrapolation below number
densities of $\sim 10$ cm-3 suggested by, e.g. Crutcher et al. (2010).

Figure 2: Median magnetic field strength vs. gas number density over time.
Black points are measured values from clouds inside the Milky Way (Crutcher et
al., 2010). Amplification is visible in diffuse gas in galaxies with weaker
initial fields.

###  3.2 Overall structure

Figure 3: Visualisations of each galaxy 1 Gyr. The top row shows surface
density projections, and the bottom row shows slices of magnetic field
strength in the midplane. The Hydro case is black because there is zero
magnetic fields everywhere.

We now proceed with an examination of the state of the galaxies after 1 Gyr of
evolution. Figure 3 shows face-on visualisations of each galaxy. The top row
shows the surface densities of gas in each galaxy. The excess of gas in the
medium and strong galaxies is clearly visible. Those galaxies also appear less
fragmented, with more flocculent spiral arms and fewer superbubble holes. The
bottom row shows slices of the magnetic field strength in the midplane of the
galactic disk. Although there is some amplification of the fields during the
galaxies evolution, the fields do not saturate at the same level and the case
with the stronger initial fields still has the strongest magnetic fields at
the end of the run. The stronger magnetic fields have less structure in them,
mostly due to them having lower star formation rates (see section 3.3), but
also due to the weaker fields being less dynamically important and having less
resistance to being pushed around by motions of the gas. Large scale spiral
structure is reflected in the magnetic fields, and the stronger field cases
result in spiral morphology that is less disrupted by superbubbles. To
understand this better, we need to examine the distribution of star formation.

Figure 4: Star formation rate vs. time for each galaxy, excluding the central
2 kpc. Star formation rate is calculated by summing the mass of all stars
formed within each time bin and dividing by the size of the bin. For the MHD
Weak and Medium galaxies, field strengths decline over time due to increasing
field strengths. Figure 5: Top row: gas surface density as a function of
galactocentric radius at 1 Gyr. Middle row: star formation rate surface
density as a function of radius. Bottom row: star formation rate surface
density vs. gas surface density, otherwise known as the Kennicutt-Schmidt
relation. Open circles indicate points within 2 kpc of the center. Diagonal
lines indicate constant gas depletion times of 108, 109, and 1010 years, as
done in (Bigiel et al., 2008)

###  3.3 Star Formation

Figure 4 summarizes the star formation history of each galaxy. To ensure a
fair comparison between the galaxies we restrict this analysis to gas outside
of the central 2 kpc, which minimizes the differences arising from the large
difference in gas depletion in the center of the galaxies in different cases.
The Hydro and MHD Weak galaxy undergo starbursts that have SFRs reaching up to
20 M⊙ yr-1. After the initial burst, the hydro galaxy settles into a constant
SFR of $\sim 2M_{\odot}yr^{-1}$ that reduces slightly by the end due to a
decreasing gas content. The MHD Weak galaxy has an an initially elevated SFR
of 5 M⊙ yr-1, which continually decreases and reaches 1 M⊙ yr-1 by the end of
the run, lower than the hydro galaxy. The MHD Medium does initially undergo a
slight starburst, but it is delayed until $\sim$75 Myr and much smaller. Once
it does begin to form stars it quickly reaches a SFR of 4 M⊙ yr-1 before
decreasing even more quickly than the MHD Weak galaxy, and ends up with a SFR
of less than 1 M⊙ yr-1. The decreasing SFR in both the MHD Weak and Medium
galaxies is due to the amplification of the magnetic fields, As they become
stronger they become more dynamically important, they limit star formation
more effectively. The MHD Strong galaxy has very limited star formation,
remaining around 1 M⊙ yr-1 or less for its entire evolution.

The differences in the remaining gas content can be seen in the top row of
figure 5, which plots the surface density of gas versus galactocentric radius.
The stronger the fields in the galaxy, the more gas remains due to the
different star formation history. The differences are most prominent in the
galactic center but exist out to around 10 kpc, beyond which all 4 galaxies
have similar gas surface densities. The second row of figure 5 shows star
formation surface density as a function of radius, averaged over the last 100
Myr of the simulation. The stronger field galaxies have enhanced star
formation in the center 2 kpc due to having more gas remaining at this point
in time. In the outer regions, star formation only occurs out to a limited
distance. It is truncated at 14 kpc in the Hydro and MHD Weak cases, 10 kpc in
the MHD Medium case, and at 8 kpc in the MHD Strong case. The truncation of
star formation in the stronger field cases happens where they have higher
surface densities which are typically association with higher star formation
rates, indicating that the magnetic fields can completely shut down star
formation if strong enough.

The final row in figure 5 combines above data to make the well-known
Kennicutt-Schmidt plot. The surface densities were measured at 1 kpc
resolution, and then binned by radius. Star formation rate surface density is
calculated using stars that formed within the last 100 Myr in each pixel to
make a fair comparison to observational tracers of star formation rate. The
dashed lines on the plot show lines of constant consumption times of 108, 109,
and 1010 years from top to bottom. Bigiel et al. (2008) label these using
efficiencies per 108 yr of 100%, 10% and 1%, respectively. The Hydro and MHD
Weak galaxies have most of their gas at surface densities of 1-10 M⊙/pc2. They
approach a consumption time of around $2\times 10^{9}$ years at 10 M⊙/pc2 with
a reduced star formation at lower surface densities. This behaviour is typical
of the Bigiel et al. (2008) data. The stronger magnetic fields push the turn
down further right, with star formation being drastically reduced at low
surface densities. They extend to higher surface densities towards the center
of the galaxies. We note that the open points represent points with radii less
than 2 kpc. The MHD Medium and Strong galaxies both have magnetic fields that
are strong enough to support the gas with levels of star formation that are
below the observations. This result is reinforced when we explicitly examine
gas support in the next section.

###  3.4 Gas Properties and Distribution

The Hydro and MHD Weak galaxies end up with higher star formation rates than
the other two galaxies despite having less gas content remaining at the end.
Our star formation model depends on the amount of dense gas. Thus we expect to
see less dense gas when there is reduced star formation. Figure 6 shows a
histogram of the total gas mass vs. number density excluding the central 2 kpc
at the final snapshot. We see progressively less star forming gas (above the
number density of 100) in the stronger field cases as expected. The left side
of the distribution is also more extended with the weaker fields. This is
because of the stronger feedback and larger bubbles. The regular small bumps
in the plot are a result of the refinement strategy of Ramses.

Figure 6: Histogram of gas number densities excluding central 2 kpc. Each bin
represents the total mass at that density inside of a disk with radius 15 kpc,
and 1 kpc high, normalized by the total gas mass of that disk. Magnetic fields
limit the amount of star forming gas that is created by narrowing the
distribution. Figure 7: Example of a typical supernova bubble forming in the
simulation. Color shows the number density of gas, with the magnetic field
lines visualized by a line-integral convolution. Overplotted are the velocity
of the gas (white), and newly formed star particles (blue). This particular
bubble occurs at a radius of 7 kpc, at t=500 Myr in the MHD Medium galaxy. As
the bubble expands, magnetic fields lines are dragged with the gas, and resist
expansion.

As individual superbubbles form due to supernova feedback, the magnetic fields
lines are dragged by the gas in the explosions, wrapping themselves around the
bubbles. Figure 7 shows a typical example of such an event. This particular
bubble has stopped expanding and does not end up escaping out of the disk.
These explosions are a major source of turbulence in the gas and play a major
role in the evolution of both the galaxy and its magnetic fields. Due to the
magnetic fields resisting being bent and compressed, they will counteract the
expansion of the bubbles as the fields are dragged. Thus the field strength
also affects the visual morphology by limiting the number and size of holes in
the gas distribution as seen in figure 3.

The gas volume distribution is affected more dramatically affected by the
magnetic fields than the mass distribution. In figure 8, we plot the volume
fraction of each phase of gas over time. We define the three phases as cold
gas (T < 5000 K), warm gas (5000 K < T < 50000 K), and hot gas (T > 50000 K).
To first order, the volume fraction is explained by the star formation; high
star formation rates lead to more supernovae which create more hot gas in
superbubbles. As the star formation decreases in the MHD Weak and Medium
galaxies, the volume fraction of hot gas decreases correspondingly. But that
is not the whole picture; the cases with magnetic fields have a lower hot
volume fraction for a given star formation rate. Between 500 -700 Myr, the MHD
Weak galaxy has roughly the same star formation rate as the Hydro galaxy, but
a systematically lower volume fraction. Similarly, when the MHD Medium galaxy
peaks in star formation around t=150 Myr, which is higher than the hydro
galaxy, it only achieves a maximum hot volume fraction of just over 20
percent. The stars that do form in the MHD Strong galaxy are hardly able to
make any bubbles at all, with the entire galaxy’s volume being dominated by
warm phase gas. The MHD Medium galaxy ends up in a similar state by the end of
the simulation.

This strongly suggests that the magnetic fields are limiting the growth of
superbubbles, as seen in figure 7. Another possible cause of the difference in
bubble volume is the clustering of stars, if the star formation is more
clustered, the supernovae combine more efficiently and will grow even larger
(Nath & Shchekinov, 2013; Keller et al., 2014). We have confirmed that the
masses of the star clusters do not change between the four galaxies, all
having similar distributions. It is commonly estimated that hot gas occupies
$\sim 50$ % of the ISM by volume (Tielens, 2010). Only the MHD Weak and Hydro
cases reflect this. The MHD Medium and Strong filling factors in figure 8
extremely low at late times.

Figure 8: Volume fraction of each phase of gas in the ISM over time. Volume is
calculated from a disk of radius 15 kpc, and 1 kpc high. Majority of the
volume is split between the warm and hot phase gas. With increasing field
strength the volume fraction of the warm phase gas increases, and hot phase
gas decreases.

Figure 9 shows the mass weighted average of the height of the gas vs. radius.
All of the galaxies have a value less than 100 pc in the center, increasing as
the disk flares outwards. The MHD Weak Galaxy has the thinnest disk, due to
its reduced star formation rate. It’s thickness is fairly similar to the Hydro
galaxy, and they both reach 250 pc by a radius of 15 kpc. The MHD Medium and
Strong galaxies are both systematically thicker, despite their drastically
reduced star formation rates. The MHD Strong galaxy is the thickest, reaching
a height of 400 pc. Because of the mass weighting this measure of thickness is
mostly set by the high density gas, but in the galaxies with stronger fields
there is less cold gas so the increased height is largely due to warm diffuse
gas which is magnetically supported. To quantify these results, we need to
examine the different contributions to the supporting pressure.

Figure 9: Mass weighted average |z| (which represents the thickness of the
disk), as a function of galactocentric radius. In each radial bin, gas within
4 kpc vertically of the midplane was included. The MHD Weak galaxy has a
thinner disk than the hydro case, but the more strongly magnetized galaxies
are thicker.

The differences in gas properties and star formation rate ultimately arise due
to the different forces (or pressures) acting on the gas, including thermal
pressure, magnetic forces, turbulent motions and gravity (Benincasa et al.,
2016). Galaxies that have stronger magnetic fields have reduced turbulence due
to the reduced stellar feedback.

Figure 10 shows mass-weighted pressures providing vertical support within 500
pc of the disk midplane versus radius, averaged over the final 100 Myr. The
curves shows thermal pressure, $P_{\textrm{thermal}}=n\,k\,T$, magnetic
pressure $P_{\textrm{B}}=B^{2}/8\pi$, and a turbulent pressure, estimated
using $P_{\textrm{turb}}=\rho v_{z}^{2}$. The top row is cold gas, the middle
row is warm gas, and the bottom row is hot gas, using the same definitions as
above (i.e. in figure 8). In the cold gas, the MHD Weak, Medium, and Hydro
galaxies have roughly the same total pressure, however with increasing
magnetic field strength there is a smaller contribution from turbulence. The
central regions of the MHD Medium and Strong galaxies have higher support to
hold up the extra gas remaining there, and the radial extent of the cold phase
is reduced. We draw attention to the difference between the Hydro and MHD Weak
cases. The turbulent support is roughly half for MHD Weak case (with the
difference being made up for by magnetic support), mirroring the halving of
the star formation rate at this time relative to Hydro.

In the warm gas, the difference in magnetic pressure is much larger between
the MHD Weak and Medium. The increase in magnetic pressure in warm gas is
largely responsible for the increased scale height in the Medium and Strong
Galaxies. All pressures are in a rough equipartition in the MHD Weak galaxy,
but the warm phase becomes dominated by magnetic pressure in the Medium and
Strong galaxies.

The hot gas is localized in high pressure superbubbles which are dominated by
thermal pressure in all galaxies. The magnetic fields expand with the hot gas,
resulting in the lower magnetic pressure in the hot gas. There is also a high
turbulent pressure, which is likely due to high-velocity flows as opposed to
small scale turbulence.

Figure 10: Comparison of different pressures in each galaxy. Included
pressures are thermal (red), magnetic B2/8$\pi$ (blue), and turbulent $\rho
v_{z}^{2}$ (yellow). Each pressure is calculated by taking a mass weighted
average of gas within 1 kpc vertically in each radial bin. The pressures are
calculated for three separate phases, using only the gas that falls within
each temperature range: cold (T < 5000 K), warm (5000 K < T < 50000 K), and
hot (T > 50000 K).

###  3.5 Disk stability

We have shown that the formation of a cold phase and subsequent star formation
is much more limited in the stronger field cases. Gravitational instability
drives the formation of the cold phase. We can quantify this using the Toomre
Q parameter for the gas,

$Q=\frac{\kappa\sqrt{c_{s}^{2}+v_{a}^{2}+\sigma_{v}^{2}}}{\pi G\Sigma},$ (4)

which accounts for the thermal, magnetic and turbulent support discussed in
the previous section. Here $c_{s}$ is sound speed, $v_{a}$ is the Alfvén
velocity, $\sigma_{v}$ is the velocity dispersion of the gas, and $\Sigma$ is
the surface density of the gas. The literature contains several extended
versions of the Toomre Q parameter (Romeo & Falstad, 2013; Körtgen et al.,
2019; Nipoti, 2023), which include adjustments for the 3D structure of the
disk and multiple components like the stellar disk. We note that the Agora
stellar disk has quite a high velocity dispersion and thus the stellar
component is quite stable and does not contribute much to the effective Q
locally. The regions of low gas Q closely match up with the locations where
stars form, confirming that this choice of Q is a reasonable approximation.
Figure 11 shows the Toomre Q parameter for each galaxy at 1 Gyr, along with
each of the individual support terms. The Alfvén velocity and sound speed were
calculated by taking a mass weighted average of gas within 250 pc of the
midplane, and the velocity dispersion in the midplane is estimated by summing
over the differences between neighbouring cells.

The hydro galaxy again has highest turbulent support which is localized around
star forming regions. Conversely, velocity dispersion provides the least
support in the magnetic galaxies. Those same regions also contain high
temperature gas which yields high sound speeds. The prevalence of hot bubbles
drops rapidly as the field strength increase from left to right in the third
row of the figure. The thermal support is lowest in dense gas seen along
spiral arms, where the magnetic support is highest.

The fourth row of figure 11 shows the Toomre Q parameter, with red regions
indicating Q$<1$. Much of the material with lower Q values has already
collapsed so these values below 1 are not informative. Each galaxy is unstable
only in regions of high surface density along spiral arms. The galaxies with
stronger magnetic fields have dramatically fewer regions that are unstable and
at large radii become totally stable against collapse, explaining the
locations of star formation in the fifth row. Our MHD Weak galaxy has
similarly placed but smoother unstable regions when compared to the hydro
galaxy. This is directly related to the addition of magnetic support shown in
the top panel. The spiral features are thus smoother and more continuous,
which is a general feature in the magnetized cases.

As discussed in Körtgen et al. (2019), magnetized galaxies may also be able to
fragment via the Parker instability rather than gravitational instability,
potentially allowing collapse in regions where $Q>1$. This would be difficult
to observe in our galaxies because stellar feedback is constantly stirring the
gas and may disrupt the wavelike structure required. However, we see
relatively little star formation outside of Toomre unstable regions. In
addition, the timescale for these large-scale Parker instabilities tends to be
long compared to turbulent crossing times.

In the outer regions of the disks, much of the gas is marginally stable
($1\leq Q\leq 2$). In this regime, the development of spirals is still
possible due to swing amplification (Binney & Tremaine, 2008). These spirals
form but do not necessarily fragment into clouds. In the MHD strong case no
dense structure forms in the outer regions (see figure 3) due to the enhanced
magnetic support, and the star formation is heavily centrally concentrated, as
shown in the bottom row.

Figure 11: Toomre Q parameter and the three support terms included in it as
shown in Equation 4. Top row: Alfvén velocityy. Second row, velocity
dispersion. Third row: sound speed. Fourth row: Toomre Q parameter. Red
regions show Q<1, which are unstable to collapse. Q values of 1-2 are unstable
to forming spirals. Fifth row: Surface density of star formation using stars
formed within the last 100 Myr.

##  4 Discussion

We simulated four galaxies with varying magnetic field strengths in order to
understand the impact on their ISM and star formation regulation. We saw clear
differences between the galaxies in their star formation rates. This result
arises due to the consumption of gas fuel and the growth of the magnetic
fields strength over time. Star formation rates can be understood via the gas
support in each case. This is clearly reflected in the stability of the gas.
Stronger fields are able to support the gas without feedback, particularly at
large radii. Weak or absent fields must rely on turbulent support associated
with star formation. Every case is internally self-consistent. However, the
combination of the star formation rates, field strengths and gas morphology
heavily favour the steady state generated by the MHD Weak case as the best
match to observed disk galaxies.

A lack of star formation also means a lack of cold, dense gas precursors for
star formation, which is reflected in the ISM gas phase results. This is tied
to the Schmidt-law star formation model, that always forms stars at a fixed
efficiency per free-fall time once gas is above a density threshold. While
this is a standard model choice, it is important to recognize that the actual
efficiency of star formation in dense gas is the subject of vigorous
theoretical and observational study. In particular, the presence of unresolved
magnetic fields and turbulence could affect this efficiency in real galaxies,
whereas the efficiency is fixed here.

It is difficult to compare simulated field strengths to observed values. We
argue that volume weighted fields are similar to observed synchrotron
estimates. In particular, the MHD Weak case creates a flat radial profile.
While this is compelling relative to the other cases, mock observations would
provide more detailed insights (Ponnada et al., 2022, 2023).

The MHD Weak galaxy was the only one that reproduced the steep power law
increase in field strength at high densities inferred from Zeeman observations
(Crutcher et al., 2010). Our ideal MHD simulated field strengths should be
upper limits, as turbulent ambipolar diffusion will reduce the field at high
densities (Heitsch et al., 2004). Our simulations also show no evidence for
the commonly inferred constant field strength in lower density gas. This
behaviour has been observed by several independent simulation groups (Körtgen
et al., 2019; Rieder & Teyssier, 2017b; Ponnada et al., 2022). Thus there is a
tension between observations and theory on this point.

It is tempting to ask which of the four galaxies are the most realistic.
Magnetic fields certainly exist in galaxies, ruling out the hydro galaxy.
Similarly, the MHD Strong galaxy had initial fields above what we infer from
observational constraints. Both it and the MHD medium case failed to match
most observed properties for the field itself, gas morphologies and star
formation.

The MHD Weak case did well on all these measures. We expect that any case
starting with a sufficiently high plasma beta ($\gtrsim 100$) would evolve to
a saturated state similar to the MHD Weak case. The main difference would be
the time it takes to do so. The key requirement is that the fields are allowed
to amplify naturally to the point of saturation.

We note that the MHD Weak galaxy achieves similar star formation rates and gas
distributions to the Hydro galaxy at intermediate times. Eventually, the gas
support includes a magnetic component similar to the turbulence and the star
formation rate dips below the hydro case. In future work, we will explore the
morphology of the field and how it is linked to its ability to support the gas
in place of stellar feedback. This may be related to the split between
turbulent and mean fields. It would also be worth testing different
morphologies in the initial condition. We started with a purely toroidal field
which could possibly bias the final field configuration.

This work does not account for cosmic rays, which could contribute
significantly to gas support in principle (e.g. Semenov et al., 2021). We note
that there is a lot of uncertainty regarding cosmic rays coupling to the gas
and how best to model them numerically. While the current work uses simple
approaches to isolate key effects, cosmic rays would be an interesting future
direction.

We also do not account for non-ideal MHD effects. Ideal MHD is a good
approximation for most of the ISM. Effects such as ambipolar diffusion may
become important in high density star forming regions. We anticipate that
these effects are comparable to the impact of the magnetic field on small star
formation generally, which is not addressed in our simple star formation
prescription. We also avoided more elaborate feedback models. Complex feedback
models would also be affected by unresolved turbulence and magnetic fields.
Though we chose simple feedback by design to keep interpretation simple, it is
worth keeping in mind that factors like star formation efficiency could be
different in different cases (e.g. with stronger magnetic fields).

Studying the properties of star forming regions and molecular cloud analogues
is a natural extension of this work to smaller scales. Magnetic fields may
affect the shapes and sizes of star forming clouds, which could also affect
small-scale star formation. To do this, we need higher resolution simulations.
We have performed zoom-in simulations, beginning with the galactic setup in
this work, seeking to characterize the affect of galactic-scale dynamics on
star forming regions (Zhao et al 2023, in prep).

Another population of interest within galaxies is that of superbubbles. JWST
observations have made it possible to perform a detailed census of
superbubbles in nearby galaxies (Watkins et al., 2023). Our results show that
magnetic fields dramatically affect the evolution of superbubbles. It would be
interesting to characterize their populations in these and higher resolution
simulations as a further constraint on how well the simulated magnetic fields
match those in real galaxies.

##  5 Conclusions

We summarize our conclusions as follows:

  * •

We have demonstrated that evolving a simulated galaxy towards a realistic
self-regulating ISM and correspondingly realistic magnetic fields is
achievable with an initially weak toroidal field and simple models for star
formation and feedback. conversely, starting with field strengths typical for
the ISM does not produce realistic end states.

  * •

Stronger magnetic fields generally reduce a galaxy’s star formation. As
magnetic fields amplify over time and their strengths increase, star formation
rates will decrease. The field strengths are also linked to the star formation
and feedback and provide an additional form of self-regulation in galaxies.

  * •

In our isolated galaxies, dynamo amplification mainly occurs in the warm
diffuse medium of the ISM, from number densities of 0.1 to 10 cm-3. The
galactic magnetic fields first saturate at high densities but continue to grow
at intermediate densities over longer timescales. Starting with high initial
fields may lock the galaxy into unrealistic field configurations.

  * •

Stronger magnetic fields generally result in reduced turbulence. This is due
to both the reduced supernovae feedback, but also from magnetic fields
limiting superbubble growth.

  * •

Magnetic fields play an important role in the vertical pressure support in
galaxies, reaching equipartition levels with turbulence in the cold phase, and
with both thermal pressure and turbulence in the warm phase. This adds an
extra component the established pressure regulation picture (Ostriker et al.,
2010; Benincasa et al., 2016).

  * •

Magnetic fields change the distribution of unstable regions in the disk, which
are well characterized with a modified gas Toomre Q. Strongly magnetized
galaxies can be completely stabilized against collapse for much of the disk.

  * •

Strong magnetic fields can dramatically reduce the size of supernova bubbles.
This can be seen visually in the simulated galaxies, and by the reduced volume
fraction of hot gas for a given star formation rate.

## Acknowledgements

The authors would like to thank Ralph Pudritz for many useful dicussions.
Hector Robinson is supported by an NSERC postgraduate scholarship, and James
Wadsley is supported by Discovery Grants from NSERC of Canada. Computational
resources for this project were enabled by a grant to James Wadsley from
Compute Canada/Digital Alliance Canada and carried out on the Niagara
Supercomputer.

## Data Availability

The data used this article will be shared upon reasonable request to the
corresponding author.

## References

  * Agertz et al. (2011) Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391
  * Basu & Roy (2013) Basu A., Roy S., 2013, MNRAS, 433, 1675
  * Beck & Wielebinski (2013) Beck R., Wielebinski R., 2013, in Oswalt T. D., Gilmore G., eds, , Vol. 5, Planets, Stars and Stellar Systems. Volume 5: Galactic Structure and Stellar Populations. p. 641, doi:10.1007/978-94-007-5612-0_13
  * Beck et al. (2019) Beck R., Chamandy L., Elson E., Blackman E. G., 2019, Galaxies, 8, 4
  * Benincasa et al. (2016) Benincasa S. M., Wadsley J., Couchman H. M. P., Keller B. W., 2016, MNRAS, 462, 3053
  * Benincasa et al. (2020) Benincasa S. M., Wadsley J. W., Couchman H. M. P., Pettitt A. R., Keller B. W., Woods R. M., Grond J. J., 2020, MNRAS, 499, 2028
  * Bigiel et al. (2008) Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1
  * Burkhart (2018) Burkhart B., 2018, ApJ, 863, 118
  * Cao & Li (2023) Cao Z., Li H.-b., 2023, ApJ, 946, L46
  * Chyży & Buta (2008) Chyży K. T., Buta R. J., 2008, ApJ, 677, L17
  * Crutcher et al. (2010) Crutcher R. M., Wandelt B., Heiles C., Falgarone E., Troland T. H., 2010, ApJ, 725, 466
  * Evans & Hawley (1988) Evans C. R., Hawley J. F., 1988, ApJ, 332, 659
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath et al. (2011) Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011, Phys. Rev. Lett., 107, 114504
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Fletcher et al. (2011) Fletcher A., Beck R., Shukurov A., Berkhuijsen E. M., Horellou C., 2011, MNRAS, 412, 2396
  * Geach et al. (2023) Geach J. E., Lopez-Rodriguez E., Doherty M. J., Chen J., Ivison R. J., Bendo G. J., Dye S., Coppin K. E. K., 2023, arXiv e-prints, p. arXiv:2309.02034
  * Girichidis et al. (2018) Girichidis P., Seifried D., Naab T., Peters T., Walch S., Wünsch R., Glover S. C. O., Klessen R. S., 2018, MNRAS, 480, 3511
  * Goldsmith et al. (2008) Goldsmith P. F., Heyer M., Narayanan G., Snell R., Li D., Brunt C., 2008, ApJ, 680, 428
  * Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
  * Han (2017) Han J. L., 2017, ARA&A, 55, 111
  * Heitsch et al. (2004) Heitsch F., Zweibel E. G., Slyz A. D., Devriendt J. E. G., 2004, ApJ, 603, 165
  * Hockney & Eastwood (1981) Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, MNRAS, 480, 800
  * Kainulainen et al. (2009) Kainulainen J., Beuther H., Henning T., Plume R., 2009, A&A, 508, L35
  * Keller et al. (2014) Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013
  * Kim et al. (2016) Kim J.-h., et al., 2016, ApJ, 833, 202
  * Kim et al. (2023) Kim C.-G., Kim J.-G., Gong M., Ostriker E. C., 2023, ApJ, 946, 3
  * Körtgen et al. (2019) Körtgen B., Banerjee R., Pudritz R. E., Schmidt W., 2019, MNRAS, 489, 5004
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Lopez-Rodriguez et al. (2023) Lopez-Rodriguez E., et al., 2023, ApJ, 942, L13
  * Martin-Alvarez et al. (2018) Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479, 3343
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, in AGU Fall Meeting Abstracts. pp SM51B–1295
  * Nath & Shchekinov (2013) Nath B. B., Shchekinov Y., 2013, ApJ, 777, L12
  * Nipoti (2023) Nipoti C., 2023, MNRAS, 518, 5154
  * Ostriker et al. (2010) Ostriker E. C., McKee C. F., Leroy A. K., 2010, ApJ, 721, 975
  * Padoan et al. (2012) Padoan P., Haugbølle T., Nordlund Å., 2012, ApJ, 759, L27
  * Pakmor et al. (2020) Pakmor R., et al., 2020, MNRAS, 498, 3125
  * Pattle et al. (2023) Pattle K., Fissel L., Tahani M., Liu T., Ntormousi E., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Astronomical Society of the Pacific Conference Series. p. 193
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 586, A138
  * Ponnada et al. (2022) Ponnada S. B., et al., 2022, MNRAS, 516, 4417
  * Ponnada et al. (2023) Ponnada S. B., et al., 2023, arXiv e-prints, p. arXiv:2309.04526
  * Rieder & Teyssier (2016) Rieder M., Teyssier R., 2016, MNRAS, 457, 1722
  * Rieder & Teyssier (2017a) Rieder M., Teyssier R., 2017a, MNRAS, 471, 2674
  * Rieder & Teyssier (2017b) Rieder M., Teyssier R., 2017b, MNRAS, 472, 4368
  * Riols & Latter (2019) Riols A., Latter H., 2019, MNRAS, 482, 3989
  * Robinson (2021) Robinson H., 2021, M.Sc Thesis, McMaster University, http://hdl.handle.net/11375/26150
  * Romeo & Falstad (2013) Romeo A. B., Falstad N., 2013, MNRAS, 433, 1389
  * Semenov et al. (2021) Semenov V. A., Kravtsov A. V., Caprioli D., 2021, ApJ, 910, 126
  * Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
  * Steinwandel et al. (2019) Steinwandel U. P., Beck M. C., Arth A., Dolag K., Moster B. P., Nielaba P., 2019, MNRAS, 483, 1008
  * Steinwandel et al. (2020) Steinwandel U. P., Dolag K., Lesch H., Moster B. P., Burkert A., Prieto A., 2020, MNRAS, 494, 4393
  * Steinwandel et al. (2022) Steinwandel U. P., Böss L. M., Dolag K., Lesch H., 2022, ApJ, 933, 131
  * Su et al. (2018) Su K.-Y., Hayward C. C., Hopkins P. F., Quataert E., Faucher-Giguère C.-A., Kereš D., 2018, MNRAS, 473, L111
  * Sugitani et al. (2011) Sugitani K., et al., 2011, ApJ, 734, 63
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Tielens (2010) Tielens A. G. G. M., 2010, The Physics and Chemistry of the Interstellar Medium
  * Watkins et al. (2023) Watkins E. J., et al., 2023, ApJ, 944, L24
  * Wissing & Shen (2023) Wissing R., Shen S., 2023, A&A, 673, A47
  * Wolfire et al. (2003) Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278
  * Zweibel (2017) Zweibel E. G., 2017, Physics of Plasmas, 24, 055402


# Multichromatic quantum superpositions in entangled two-photon absorption
spectroscopy

M Wittkop,a Juan M. Marmolejo-Tejada,a Martín A. Mosqueraa,∗
Department of Chemistry and Biochemistry
103 Chemistry and Biochemistry Building, Bozeman, MT-59717, USA
<EMAIL_ADDRESS>
###### Abstract

Quantum information science is driving progress in a vast number of scientific
and technological areas that cover molecular spectroscopy and matter-light
interactions in general. In these fields, the ability to generate quantum
mechanically-entangled photons is opening avenues to explore the interaction
of molecules with quantum light. This work considers an alternative way of
correlating photons by including energy superpositions. We study how the
multichromatic quantum superposition, or color superposition of photon-pair
states, influences the optical properties of organic chromophores. This work
uses electronic structure calculations based on time-dependent density
functional theory, and a simple modification of the standard entangled two-
photon absorption theory. Our calculations show that it is possible to
substantially modify the optical absorption cross section of molecules, where
constructive and destructive interferences are computed. The quantum
interference effects are more pronounced than the constructive ones. These
quantum effects, or related ones, could be observed in quantum spectroscopic
experiments where qudit photon states are generated.

##  1 Introduction

Quantum spectroscopy offers tools to elucidate molecular systems and materials
that both expand and complement techniques based on classical light [1, 2, 3].
This has motivated work that is constantly demonstrating the significant
transformative potential of quantum light to understand molecular function and
open technological opportunities. Such advances are taking place in parallel
with scientific and engineering fields such as photonic quantum computing [4,
5], where the precise and accurate control of correlated photons could bring
substantial advantages for diverse, cutting-edge applications.

A specific phenomenon that has received widespread attention recently is the
absorption and emission of entangled photon pairs by molecules. Entangled two-
photon absorption (ETPA) [6], for instance, introduces physical correlations
that are not possible in classical two-photon absorption (TPA) spectroscopy
[7]. ETPA commonly relies on the generation of entangled photon-pairs, which
usually occurs through the well-known spontaneous parametric down conversion
(SPDC) process [8]. SPDC provides photon pairs where the polarization of the
photons are quantum mechanically correlated, and such entanglement can be
verified in EPR-like (Einstein-Podolsky-Rosen) devices [9, 10, 11]. These
photons are emitted within an quantum area, and the photons in each entangled
pair are delayed with respect to one another [12, 13]. Such delay is expressed
in terms of an entanglement time, which endow ETPA techniques with unique
properties that to-date continue to be explored by the community. Entangled
photons can be emitted through other mechanisms such as molecular pathways
[14, 15], quantum dots [16, 17], or semiconducting devices [18, 19].

A family of chromophores have been investigated recently: molecules such as
Rhodamine 6G [20], Zinc-TPP [21], and thiophene dendrimers [22], among other
dyes [23]. ETPA absorption of a molecular unit is quantified commonly through
ETPA cross-sections, which have the same units as classical one-photon cross
sections, usually expressed in $\mathrm{cm}^{2}$ units. Rigorous experimental
efforts suggest that ETPA offers considerable quantum advantage over classical
TPA (CTPA) [23], especially at low photon (quantum light) fluxes (at extremely
low fluxes ETPA will dominate significantly over CTPA). However, the
estimation of entangled cross sections with very high accuracy is the subject
of current efforts [24, 23]. This is a strong motivation, in our opinion, to
further the understanding of the interaction of molecules and materials with
quantum light and their connection to quantum technologies. One can
hypothesize that experimental and theoretical techniques will continue to
advance in these directions, unlocking unexpected and highly beneficial
quantum phenomena in a wide variety of systems.

Experimental and theoretical studies so far have focused on the interplay
between single entangled photon pairs and molecules, where the frequencies of
the photons are assumed to be given. For example, in degenerate pumping the
frequency of both photons are the same. Recently, however, energy
superposition of photons has been achieved for individual photons [25, 26, 27]
and for photon pairs [18]. Energy superpositions give rise to the well-known
“qudit” states, which generalize the concept of the qubit. Also known as color
superposition, in this phenomenon, the color of the photons is undetermined,
and each single photon is in a superposition of two (or more) colors. Even
though the interaction between photon qubit (or qudit) states and molecules
has not been reported experimentally so far, their effects can be explored
theoretically. We do so in this case for three chromophores of interest:
flavin mononucleotide (FMN) [28, 29], topotecan (TPT) [30, 31, 32, 33], and
lucifer yellow (LY) [34]. This work studies the molecular absorption of
entangled photon-pairs that also feature multichromatic superpositions and
polarization entanglement. These photon pairs could form qudits of six-fold
dimensionality (or eight-fold if four colors are used); this work, however,
focuses on a qubit representation, as specified herein, but a higher-
dimensionality of quantum states is possible. We find that the cross sections
in this case show signatures of constructive and destructive interferences,
depending on the location of the quantum superposition in the Bloch sphere.
For a special set of angles in the Bloch sphere, we notice that the
destructive interference can be quite substantial, and for other angles, we
observed that the absorption cross section can be enhanced by close to an
order of magnitude, whereas the quantum interference can lower absorption by
up to two orders of magnitude. These findings then suggest the phenomenon of
color-superposition could be of interest for additional quantum control of
ETPA cross sections and related experiments.

##  2 Theory

This work focuses on the interaction of photon pairs with vertical electronic
excited states only. This is a common approximation that is employed in CTPA
spectroscopy because full vibronic transitions are quite demanding,
computationally.

For photon-pairs characterized by two unique frequencies $\omega_{1}$ and
$\omega_{2}$, entanglement time $T_{\mathrm{e}}$, and entanglement area
$A_{\mathrm{e}}$, the absorption cross-section is given by [35]:
$\sigma_{f,0}=(4\pi^{3}\alpha a_{0}^{5})/(A_{\mathrm{e}}\tau_{0}c)\times
g(\omega_{\mathrm{T}}-\Omega_{f})|W_{f,0}|^{2}$, where $g$ is the line shape
function, $\omega_{\mathrm{T}}$ is the sum of the two photon frequencies,
$\tau_{0}$ is the atomic unit of time ($\tau_{0}=m_{e}a_{0}^{2}/\hbar$),
$\Omega_{f}$ is the excitation energy for transition from the ground state to
the excited state labeled $f$ (the ground state is labeled as the $0$-th
state), $a_{0}$ is the Bohr length, $\alpha$ the fine structure constant, and
$c$ the speed of light. The cross-section then has the units that arise from
the term $(4\pi^{3}\alpha a_{0}^{5})/(A_{\mathrm{e}}\tau_{0}c)$. This result
for $\sigma_{f,0}$ can be derived using second perturbation theory and
assuming that the photons are modeled as uniform plane wave packets with
fronts that are spatially separated by a distance of $cT_{\mathrm{e}}$. In
terms of random photon detection times, the time difference between photon
detections is then in average the entanglement time, $T_{\mathrm{e}}$.

The function $W_{f,0}$ reads:

$W_{f,0}(\omega_{1},\omega_{2},T_{e})=\sqrt{\frac{\omega_{1}\omega_{2}}{T_{e}}}S_{f,0}$
(1)

where $\omega_{1}$ and $\omega_{2}$ are the frequencies of the first and
second incoming photons, correspondingly. In the expressed equations, all the
frequencies, dipole moments, and entanglement times are expressed in atomic
units (a.u.); this includes $\kappa$ and $\Gamma$. The term $S_{f,0}$
represents the transition function (in atomic units):

$S_{f,0}=\sum_{j}\frac{(\vec{\mu}_{fj}\cdot\vec{\epsilon}_{2})(\vec{\mu}_{j0}\cdot\vec{\epsilon}_{1})}{\Omega_{j}-\omega_{1}-\mathrm{i}\kappa}\Big{\\{}1-\exp\big{[}-\mathrm{i}(\Omega_{j}-\omega_{1}-\mathrm{i}\kappa)T_{\mathrm{e}}\big{]}\Big{\\}}+(1\leftrightarrow
2)$ (2)

where $\kappa$ represents the inverse of the lifetime of the intermediate
virtual state, $\vec{\epsilon}_{1}$ and $\vec{\epsilon}_{2}$ are the
polarizations of the first and second photon, respectively. The transition
dipole vector for a transition from ground state to excited state “$j$” is
denoted as $\vec{\mu}_{j0}$, whereas $\vec{\mu}_{fj}$ denotes the transition
dipole vector for a transition from the $j$-th excited state into the final
state $f$. The lineshape function $g$ is a described in terms of a Lorentzian
profile of the form
$g(\Delta\omega)=\pi^{-1}(\Gamma/2)/[\Delta\omega^{2}+(\Gamma/2)^{2}]$.

Figure 1: Pictorial representation of the present theoretical model: A two-
photon quantum packet is directed towards a molecule. The two-photon state is
a quantum superposition of two possible entangled states, a monochromatic
pair, labeled “MC”, and a bichromatic pair, referred to as the “BC” pair. The
molecule can then be excited through different pathways where different
intermediate states are involved in the overall two-photon excitation, as
suggested in this figure.

We now consider the two-photon packet as being in a superposition of two
states: “MC”, in which the photons are monochromatic
($\omega_{1}=\omega_{2}=\omega_{\mathrm{T}}/2$), and “BC”, in which the
photons are bichromatic, $\omega_{1}^{\prime}\neq\omega_{2}^{\prime}$; the
“primed” quantities are assigned to the BC state. It is assumed that MC and BC
photon-pair quantum states have entanglement times $T_{\mathrm{e}}$ and
$T_{\mathrm{e}}^{\prime}$, respectively, but have assigned the same area
$A_{\mathrm{e}}$. The quantum superimposed two-photon state is thus described
by:

$|\Psi_{\gamma}\rangle=\cos\Big{(}\frac{\theta}{2}\Big{)}|\mathrm{MC}\rangle+\sin\Big{(}\frac{\theta}{2}\Big{)}e^{\mathrm{i}\phi}|\mathrm{BC}\rangle$
(3)

we refer to this state as a “multichromatic superposition” (MCS). The MCS is
thus controlled by the Bloch sphere parameters $\theta$ and $\phi$, where
$0\leq\phi\leq 2\pi$ and $0\leq\theta\leq\pi$. Because the MCS photon
configuration obeys the standard linear superposition of states, the function
$\tilde{W}_{f,0}$ must also transform as:

$W_{f,0}^{\mathrm{QS}}=\cos\Big{(}\frac{\theta}{2}\Big{)}W_{f,0}^{\mathrm{MC}}+\sin\Big{(}\frac{\theta}{2}\Big{)}e^{\mathrm{i}\phi}W_{f,0}^{\mathrm{BC}}$
(4)

The amplitude $W_{f,0}^{\mathrm{QS}}$, as a quantum mechanical transition
element, is also described in terms of the Bloch-sphere angles $\theta$ and
$\phi$. The $W$ functions are evaluated as
$W^{\mathrm{MC}}_{f,0}=W_{f,0}(\omega_{1},\omega_{2},T_{e})$ and
$W^{\mathrm{BC}}_{f,0}=W_{f,0}(\omega_{1}^{\prime},\omega_{2}^{\prime},T_{e}^{\prime})$.
Figure 1 shows a pictorial summary of the theoretical concept explored in this
work.

The cross section for a discrete MCS transition of interest is:

$\sigma_{0\rightarrow f}^{\mathrm{QS}}=\frac{4\pi^{3}\alpha
a_{0}^{5}}{A_{e}\tau_{0}c}g(0)|W_{f,0}^{\mathrm{QS}}|^{2}$ (5)

where $\omega_{f,0}=\omega_{T}$, so $g(0)=2/\pi\Gamma$. The above equation is
a relatively simple extension of the standard formula used in standard ETPA
spectroscopy, but it now includes the possibility of there being color-
superposition. This formula supposes, as expected, that dissipation effects
and thereby decoherence in the generation of these MCS photon states are
minimal. It is important to notice that perturbation theory demands energy
conservation, as the interaction time is assumed to last for a very long time.
For this reason we have that
$\omega_{1}+\omega_{2}=\omega_{1}^{\prime}+\omega_{2}^{\prime}=\omega_{\mathrm{T}}$.

##  3 Computational Method

Figure 2: Molecular structures considered in this work: Flavin mononucleotide,
topotecan, and lucifer yellow.

The transition dipoles are calculated by means of linear-response time-
dependent density functional theory (TDDFT). We use the so-called unrelaxed
dipoles, or dipoles that come from “CIS-like” linear response TDDFT
wavefunctions. This is a practical approximation that works reasonably well,
and is computationally efficient. We extract these dipoles with an in-house
code that is based on the quantum chemistry suite NWChem [36], version 7.0. We
use the standard B3LYP exchange-correlation functional, and the 6-31G* basis
set, which is commonly applied to determine excited-state transitions in the
optical regime.

As discussed before, there is yet a source ETPA enhancement that is not well-
understood from a theory perspective. In Ref. [37], the use of radiative
lifetimes was suggested to obtain ETPA cross-section values in the range of
what has been experimentally observable. Motivated by the common practice in
TPA theory of using a standard value for the broadening, for the ETPA
calculations, we assume a value of $\Gamma$ that corresponds to
$10^{-8}~{}\mathrm{eV}$, reflecting a final state lifetime of approximately
$1.0~{}\mathrm{ns}$. For the intermediate state we assume
$\kappa=0.01~{}\mathrm{eV}$. In all the ETPA calculations the entanglement
area is $A_{e}=1.0\times 10^{-8}~{}\mathrm{cm}^{2}$. For the standard ETPA
(non-color superposition), we assume $T_{\mathrm{e}}=100~{}\mathrm{fs}$. This
same value, $T_{\mathrm{e}}=100~{}\mathrm{fs}$, is applied for the MCS (color
superimposed) ETPA, for the BC and MC quantum states of TPT and LY (Figure 2)
we take $T_{\mathrm{e}}=T_{\mathrm{e}}^{\prime}=100~{}\mathrm{fs}$, but for
FMN we set $T_{\mathrm{e}}=100~{}\mathrm{fs}$ and
$T_{\mathrm{e}}^{\prime}=75~{}\mathrm{fs}$ ($T_{\mathrm{e}}^{\prime}$ being
the entanglement time of the BC quantum state).

Classical TPA spectra are computed for comparison as well. For these we take
the standard final state broadening factor, the classical equivalent of
$\Gamma$, as $0.1~{}\mathrm{eV}$, and intermediate linewidth, the analogue of
$\kappa$, as $0.05~{}\mathrm{eV}$. We assume that the photons have cross-
polarization: so if one photon has horizontal polarization, the other has
vertical polarization. The cross sections are averaged with respect to all
possible molecular orientations [38].

##  4 Results and Discussion

Figure 3: Entangled and classical two-photon absorption cross sections: a)
classical TPA profiles of FMN and TPT; b), monochromatic ETPA cross section as
a function of $\omega_{\mathrm{h}}$ for FMN and TPT; c), ETPA cross section
vs. entanglement time for FM, TPT, and LY at energies 1.64, 1.94, and 2.33 eV,
respectively. d), ETPA profile of LY, and CTPA cross section (inset).

To illustrate the proof-of-concept for the effect of color/multichromatic
superpositions on the ETPA cross sections, we first compute the conventional
TPA spectrum and the standard ETPA spectra. Then, we proceed to examine the
MCS-ETPA properties of the chromophores selected for our theoretical study.
These chromophores are flavin mononucleotide (FMN), topotecan (TPT), and
lucifer yellow (LY); their molecular drawings are shown in Figure 2. FMN is a
biomolelcule produced from riboflavin that forms part of NADH Hydrogenase, TPT
is used as a chemotherapy drug, and LY is utilized in spectroscopic studies.
The spectra in this work are reported in terms of half the frequency of the
total two-photon frequency. This half-frequency is denoted as
$\omega_{\mathrm{h}}=\omega_{\mathrm{T}}/2$.

The theoretical CTPA spectra of FMN, TPT, and LY, are shown in Figure 3.a. and
d. (as an inset plot). We note that FMN has a classical TPA cross section of
about 220 GM at around a half-frequency of 1.64 eV (760 nm). As is commonly
seen in theoretical TPA spectra, the classical TPA shows comparable cross-
section values at higher $\omega_{\mathrm{h}}$ frequencies. A similar
numerical trend is noticed for TPT, where the CTPA cross section is about 188
GM at 1.94 eV (640 nm), and then it slightly decreases to about 41 GM at 2.12
eV, but it raises again to 70 GM around $\omega_{\mathrm{h}}=2.3$ eV. LY,
however, shows a different profile, as can be seen in Figure 3.d, where it is
relatively low - under 1 GM for half-frequencies between 2.1 and 2.7 eV, where
it raises to about 8 GM at 2.9 eV.

Figure 4: Purely bichromatic ETPA. The frequencies of BC quantum state are
$\omega_{1}^{\prime}=1/3\times\omega_{T}$,
$\omega_{2}^{\prime}=2/3\times\omega_{T}$: a), BC ETPA spectra of FMN and TPT,
and their comparison to the monochromatic counterpart; b), variation of FMN
and TPT BC ETPA cross section as a function of the entanglement time, the
frequencies $\omega_{\mathrm{h}}$ are the same as in Figure 3.c; c), BC ETPA
spectrum of LY, inset shows dependency on $T_{\mathrm{e}}$ at 2.9 eV.

Regarding standard monochromatic ETPA
($\omega_{1}=\omega_{2}=\omega_{\mathrm{h}}$), Figure 3.b shows the entangled
cross section profiles of FMN and TPT, and Figure 3.d shows that of LY. For
these three systems we use $T_{\mathrm{e}}=100~{}\mathrm{fs}$. For these
systems we see a correlation between the CTPA and ETPA spectra, however, for
TPT such connection between classical and entangled TPA is somewhat less
evident, as there are some changes between relative peak heights. FMN and TPT
have ETPA cross section values in the same order of magnitude,
$10^{-20}~{}\mathrm{cm}^{2}$. The maximum values of FMN and TPT occur at the
same frequencies as their classical TPA counterparts, with values close to
$3.0\times 10^{-20}~{}\mathrm{cm}^{2}$ and $1.8\times
10^{-20}~{}\mathrm{cm}^{2}$, respectively. As mentioned earlier, for the three
systems considered we use the same linewidth factors. The entanglement time is
another variable that affects the magnitude of the ETPA cross sections; the
effect of varying this is displayed in Figure 3.c, where we observe, as
expected, that the cross section can vary by one order of magnitude, or more,
if the $T_{\mathrm{e}}$ is further increased beyond the plot time range. For
the LY chromophore, we observe theoretical cross section values lower than
those of FMN and TPT, which indicate ETPA activity around 2.3 eV (540 nm) and
2.5 eV (496 nm), and then higher values at 2.9 eV (430 nm).

Figure 4 displays the standard ETPA spectra under a non-degenerate condition,
where the two photons have different frequencies. In this case we choose
$\omega_{1}^{\prime}=\omega_{T}/3$, and
$\omega_{2}^{\prime}=2\omega_{\mathrm{T}}/3$ and refer to this as “bichromatic
ETPA”. The bichromatic ETPA spectrum of each molecule shows a few differences
with respect to the monochromatic ones. These are mainly a slight enhancement
of the ETPA cross section, up to 30 %. However, between 2.0 and 2.4 eV, the
degenerate-pumping ETPA spectrum of TPT is not improved much by the
bichromatic condition. As in the degenerate case, the cross section at low
photon frequency varies by about one order of magnitude for $T_{\mathrm{e}}$
in the range 20 - 100 fs. Due to the parameter $\kappa$, the ETPA cross-
sections in all cases have their oscillations damped. The oscillatory pattern
of the BC ETPA case is slightly different than in the monochromatic case, but
not too significantly.

Figure 5: ETPA cross sections based on quantum color superpositions, or
quantum MCS. The Bloch sphere parameters are fixed an taken as
$\theta=60^{\circ}$ and $\phi=0^{\circ}$. For TPT and LY the entanglement
times of the BC and MC states are the same
$T_{\mathrm{e}}=T_{\mathrm{e}}^{\prime}=100~{}\mathrm{fs}$, but for FMN we
study $T_{\mathrm{e}}=100~{}\mathrm{fs}$ and
$T_{\mathrm{e}}^{\prime}=75~{}\mathrm{fs}$. a) Shows the standardard quantum
MCS ETPA spectra of FMN, b), that of TPT, and c), LY.

Having determined the difference between the separate MC and BC ETPA cases, we
now proceed to examine the multichromatic quantum superposition effect, as
expressed in Equation (4). We focus on highest intensity transitions, which
happen at 1.64, 1.94, and 2.9 eV for FMN, TPT, and LY, respectively. Figure
5.a shows the comparison between the quantum superposition effect and the
standard MC and BC ETPA cases where there is no superposition. For a polar
angle of $\theta=60^{\circ}$, the two-photon state is a quantum mixture of 75
% the monochromatic state, and 25 % the bichromatic state, yet we see that the
quantum cross section is over 100 % enhanced with respect to the monochromatic
ETPA case, for a frequency of $1.64$ eV. A similar enhancement takes place for
TPT and LY (Figures 5.b and .c, respectively), at 1.94 eV and 2.9 eV,
correspondingly. This enhancement is also due to setting the phase factor as
$\phi=0^{\circ}$. In additional preliminary calculations, we have noted
$\phi=0^{\circ}$ leads to quantum constructive effects.

For the three molecular systems, FMN, TPT, and LY, Figure 6 shows the
variation of cross-section values with respect to the Bloch sphere parameters
(again, we examine this at the frequencies 1.64 eV, 1.94 eV, and 2.9 eV,
correspondingly). The most interesting behavior in terms of quantum
constructive effects takes place at the points where theta is between
$60^{\circ}$ and $120^{\circ}$, while $\phi$ being close to $0^{\circ}$ or
$180^{\circ}$. There is thus a region of constructive interference, where the
enhancement in cross section doubles. A higher enhancement is possible for
different combinations of entanglement times [39], or if the quantum
superposition has an additional constructive effect that improves the
lifetimes of the intermediate and final excited states involved in the ETPA
quantum process (which are not explored in this work). On the other hand, a
significant destructive quantum interference takes place in a region center
around $\theta=90^{\circ}$ and $\phi=180^{\circ}$. The cross section value
drops by up to two orders of magnitude under the current parameter selection.
It drops by two orders of magnitude (from $10^{-20}~{}\mathrm{cm}^{2}$ to
$10^{-22}~{}\mathrm{cm}^{2}$) for FMN, one order of magnitude for TPT, and two
for LY. Figure 6.b shows the dependency of $\sigma$ with respect to $\theta$
for $\phi=90^{\circ}$. Clearly, the quantum superposition engenders behaviors
that are not possible under the separate circumstances, which is to be
expected from this type of phenomenon. Therefore, besides present in coherent
electronic transfer [40], destructive quantum interference is also possible in
photon absorption.

Figure 6: Two-dimensional heat plots for FMN (subfigure a) at
$\omega_{\mathrm{h}}=1.64~{}\mathrm{eV}$, TPT (subfigure c) at 1.94 eV, and LY
(subfigure d) at 2.9 eV, and a “cut” for $\phi=90^{\circ}$ for the FMN system.

As is common in quantum coherent phenomena, there are conditions that must be
satisfied to observe quantum superposition effects. On of them is the
suppression of decoherence related to vibrational degrees of freedom. These
occur during the emission of entangled photons and the absorption of these. It
is currently a subject of intensive research determining the influence of
entangled photons on the lifetime of excited states. But there are indications
that such lifetimes are enhanced by entangled photon absorption [41]. In this
work, we assumed fixed linewidth parameters in line with values used before
[37] to obtain theoretical values consistent with experimental measurements.
However, as in standard ETPA, the effect of MCS on the lifetimes could be a
subject of further research as well, as additional quantum correlations could
have unexpected consequences on electronic lifetimes and other quantum
properties.

##  5 Conclusion

In this work we investigated the interaction between entangled photon pairs,
that in addition to having the standard polarization correlation, feature
energy superpositions. That is, the frequency of the photons are undetermined
prior to interaction with matter. We referred to this phenomenon as a
multichromatic superposition. The state of the entangled photon pair was
represented in the well-known Bloch sphere, or qubit space. In comparison to
standard ETPA simulations, we observed the emergence of constructive and
destructive effects in the quantum cross section profiles. The enhancement (or
constructive) of cross section was computed to be improved by nearly a 100 %,
whereas the desctructive interference could reduce the cross section by around
two orders of magnitude. Our work then suggests that these types of coherent
effects could be added to the toolkits of quantum control within the context
of optical quantum spectroscopy.

##  6 Acknowledgments

The authors kindly thank the MonArk NSF Quantum Foundry supported by the
National Science Foundation Q-AMASE-i program under NSF award No. DMR-1906383.

## References

  * [1] Audrey Eshun, Oleg Varnavski, Juan P Villabona-Monsalve, Ryan K Burdick, and Theodore Goodson III.  Entangled photon spectroscopy.  Accounts Chem. Res., 55(7):991–1003, 2022.
  * [2] Ying-Zhong Ma and Benjamin Doughty.  Nonlinear optical microscopy with ultralow quantum light.  J. Phys. Chem. A, 125(40):8765–8776, 2021.
  * [3] Bahaa EA Saleh, Bradley M Jost, Hong-Bing Fei, and Malvin C Teich.  Entangled-photon virtual-state spectroscopy.  Phys. Rev. Lett., 80(16):3483, 1998.
  * [4] Emanuele Pelucchi, Giorgos Fagas, Igor Aharonovich, Dirk Englund, Eden Figueroa, Qihuang Gong, Hübel Hannes, Jin Liu, Chao-Yang Lu, Nobuyuki Matsuda, et al.  The potential and global outlook of integrated photonics for quantum technologies.  Nature Reviews Physics, 4(3):194–208, 2022.
  * [5] Sergei Slussarenko and Geoff J Pryde.  Photonic quantum information processing: A concise review.  Appl. Phys. Rev., 6(4):041303, 2019.
  * [6] Frank Schlawin, Konstantin E Dorfman, and Shaul Mukamel.  Entangled two-photon absorption spectroscopy.  Accounts Chem. Res., 51(9):2207–2214, 2018.
  * [7] Oleg Varnavski and Theodore Goodson III.  Two-photon fluorescence microscopy at extremely low excitation intensity: The power of quantum correlations.  J. Am. Chem. Soc., 142(30):12966–12975, 2020.
  * [8] Chao Zhang, Yun-Feng Huang, Bi-Heng Liu, Chuan-Feng Li, and Guang-Can Guo.  Spontaneous parametric down-conversion sources for multiphoton experiments.  Advanced Quantum Technologies, 4(5):2000132, 2021.
  * [9] DN Klyshko.  Combine epr and two-slit experiments: Interference of advanced waves.  Phys. Lett. A, 132(6-7):299–304, 1988.
  * [10] DN Klyshko.  A simple method of preparing pure states of an optical field, of implementing the einstein–podolsky–rosen experiment, and of demonstrating the complementarity principle.  Soviet Physics Uspekhi, 31(1):74, 1988.
  * [11] Áulide Martínez-Tapia, Samuel Corona-Aquino, Chenglong You, Rui-Bo Jin, Omar S Magaña-Loaiza, Shi-Hai Dong, Alfred B U’Ren, and Roberto de J León-Montiel.  Witnessing entangled two-photon absorption via quantum interferometry.  arXiv preprint arXiv:2208.11387, 2022.
  * [12] YH Shih and AV Sergienko.  A two-photon interference experiment using type ii optical parametric down conversion.  Phys. Lett. A, 191(3-4):201–207, 1994.
  * [13] AV Sergienko, YH Shih, and MH Rubin.  Experimental evaluation of a two-photon wave packet in type-ii parametric downconversion.  JOSA B, 12(5):859–862, 1995.
  * [14] Mohammad Rezai, Jörg Wrachtrup, and Ilja Gerhardt.  Polarization-entangled photon pairs from a single molecule.  Optica, 6(1):34–40, 2019.
  * [15] Pietro Lombardi, Maja Colautti, Rocco Duquennoy, Ghülam Murtaza, Prosenjit Majumder, and Costanza Toninelli.  Triggered emission of indistinguishable photons from an organic dye molecule.  Appl. Phys. Lett., 118(20):204002, 2021.
  * [16] Francesco Basso Basset, Francesco Salusti, Lucas Schweickert, Michele B Rota, Davide Tedeschi, Saimon F Covre da Silva, Emanuele Roccia, Val Zwiller, Klaus D Jöns, Armando Rastelli, et al.  Quantum teleportation with imperfect quantum dots.  npj Quantum Inf., 7(1):7, 2021.
  * [17] Daniel Huber, Marcus Reindl, Johannes Aberl, Armando Rastelli, and Rinaldo Trotta.  Semiconductor quantum dots as an ideal source of polarization-entangled photon pairs on-demand: a review.  J. Optics, 20(7):073002, 2018.
  * [18] Michael Kues, Christian Reimer, Piotr Roztocki, Luis Romero Cortés, Stefania Sciara, Benjamin Wetzel, Yanbing Zhang, Alfonso Cino, Sai T Chu, Brent E Little, et al.  On-chip generation of high-dimensional entangled quantum states and their coherent control.  Nature, 546(7660):622–626, 2017.
  * [19] Davide Grassani, Stefano Azzini, Marco Liscidini, Matteo Galli, Michael J Strain, Marc Sorel, JE Sipe, and Daniele Bajoni.  Micrometer-scale integrated silicon source of time-energy entangled photons.  Optica, 2(2):88–94, 2015.
  * [20] Dmitry Tabakaev, Matteo Montagnese, Geraldine Haack, Luigi Bonacina, J-P Wolf, Hugo Zbinden, and RT Thew.  Energy-time-entangled two-photon molecular absorption.  Phys. Rev. A, 103(3):033701, 2021.
  * [21] Ryan K Burdick, George C Schatz, and Theodore Goodson III.  Enhancing entangled two-photon absorption for picosecond quantum spectroscopy.  J. Am. Chem. Soc., 143(41):16930–16934, 2021.
  * [22] M R Harpham, Ö Süzer, Chang-Qi Ma, Peter Bäuerle, and Theodore Goodson III.  Thiophene dendrimers as entangled photon sensor materials.  J. Am. Chem. Soc., 131(3):973–979, 2009.
  * [23] Kristen M Parzuchowski, Alexander Mikhaylov, Michael D Mazurek, Ryan N Wilson, Daniel J Lum, Thomas Gerrits, Charles H Camp Jr, Martin J Stevens, and Ralph Jimenez.  Setting bounds on entangled two-photon absorption cross sections in common fluorophores.  Phys. Rev. Appl., 15(4):044012, 2021.
  * [24] Alexander Mikhaylov, Ryan N Wilson, Kristen M Parzuchowski, Michael D Mazurek, Charles H Camp Jr, Martin J Stevens, and Ralph Jimenez.  Hot-band absorption can mimic entangled two-photon absorption.  J. Phys. Chem. Lett., 13(6):1489–1493, 2022.
  * [25] Eva Zakka-Bajjani, François Nguyen, Minhyea Lee, Leila R Vale, Raymond W Simmonds, and José Aumentado.  Quantum superposition of a single microwave photon in two different’colour’states.  Nat. Phys., 7(8):599–603, 2011.
  * [26] Philipp Treutlein.  Photon qubit is made of two colors.  Physics, 9:135, 2016.
  * [27] Stéphane Clemmen, Alessandro Farsi, Sven Ramelow, and Alexander L Gaeta.  Ramsey interference with single photons.  Phys. Rev. Lett., 117(22):223601, 2016.
  * [28] Juan P Villabona-Monsalve, Oleg Varnavski, Bruce A Palfey, and Theodore Goodson III.  Two-photon excitation of flavins and flavoproteins with classical and quantum light.  J. Am. Chem. Soc., 140(44):14562–14566, 2018.
  * [29] Rachael J Homans, Raja U Khan, Michael B Andrews, Annemette E Kjeldsen, Louise S Natrajan, Steven Marsden, Edward A McKenzie, John M Christie, and Alex R Jones.  Two photon spectroscopy and microscopy of the fluorescent flavoprotein, ilov.  Phys. Chem. Chem. Phys., 20(25):16949–16955, 2018.
  * [30] Thomas G Burke, Henryk Malak, Ignacy Gryczynski, Zihou Mi, and Joseph R Lakowicz.  Fluorescence detection of the anticancer drug topotecan in plasma and whole blood by two-photon excitation.  Anal. Biochem., 242(2):266–270, 1996.
  * [31] Peter TC So, Chen Y Dong, Barry R Masters, and Keith M Berland.  Two-photon excitation fluorescence microscopy.  Ann. Rev. Biomed. Eng., 2(1):399–429, 2000.
  * [32] Asma Yasmeen Khan and Gopinatha Suresh Kumar.  Exploring the binding interaction of potent anticancer drug topotecan with human serum albumin: Spectroscopic, calorimetric and fibrillation study.  J. Biomol. Struct. Dyn., 36(9):2463–2473, 2018.
  * [33] Maria Rosaria di Nunzio, YiLun Wang, and Abderrazzak Douhal.  Spectroscopy and dynamics of topotecan anti-cancer drug comprised within cyclodextrins.  J Photoch. Photobio. A, 266:12–21, 2013.
  * [34] Arne S Kristoffersen, Svein R Erga, Børge Hamre, and Øyvind Frette.  Testing fluorescence lifetime standards using two-photon excitation and time-domain instrumentation: rhodamine b, coumarin 6 and lucifer yellow.  J. Fluoresc., 24:1015–1024, 2014.
  * [35] Hong-Bing Fei, Bradley M Jost, Sandu Popescu, Bahaa EA Saleh, and Malvin C Teich.  Entanglement-induced two-photon transparency.  Phys. Rev. Lett., 78(9):1679, 1997.
  * [36] Marat Valiev, Eric J Bylaska, Niranjan Govind, Karol Kowalski, Tjerk P Straatsma, Hubertus JJ Van Dam, Dunyou Wang, Jarek Nieplocha, Edoardo Apra, Theresa L Windus, and W A de Jong.  Nwchem: A comprehensive and scalable open-source solution for large scale molecular simulations.  Comput. Phys. Commun., 181(9):1477–1489, 2010.
  * [37] Gyeongwon Kang, Kobra Nasiri Avanaki, Martín A Mosquera, Ryan K Burdick, Juan P Villabona-Monsalve, Theodore Goodson III, and George C Schatz.  Efficient modeling of organic chromophores for entangled two-photon absorption.  J. Am. Chem. Soc., 142(23):10446–10458, 2020.
  * [38] PR Monson and WM McClain.  Polarization dependence of the two-photon absorption of tumbling molecules with application to liquid 1-chloronaphthalene and benzene.  J. Chem. Phys., 53(1):29–37, 1970.
  * [39] Sajal Kumar Giri and George C Schatz.  Manipulating two-photon absorption of molecules through efficient optimization of entangled light.  J. Phys. Chem. Lett., 13(43):10140–10146, 2022.
  * [40] Justin P Bergfield, Henry M Heitzer, Colin Van Dyck, Tobin J Marks, and Mark A Ratner.  Harnessing quantum interference in molecular dielectric materials.  ACS Nano, 9(6):6412–6418, 2015.
  * [41] Federica Ricci, Haraprasad Mandal, Muaaz Wajahath, Ryan K Burdick, Juan P Villabona-Monsalve, Saber Hussain, and Theodore Goodson III.  Investigations of coherence in perovskite quantum dots with classical and quantum light.  J. Phys. Chem. C, 2023.