Neutron star pulse profile observations as extreme gravity probes
Hector O. Silva Nicolás Yunes




The x-ray emission of hot spots on the surface of neutron stars is the
prime target of the Neutron star Interior Composition Explorer (NICER).

These x-ray pulse profiles not only encode information of the bulk properties of these stars,
which teaches us about matter at supranuclear densities, but also about the spacetime
curvature around them which teaches us about relativistic gravity.

We explore the possibility of performing strong-gravity tests with NICER observations
using a recently developed pulse profile model beyond general relativity.

Our results suggest that NICER can in principle place constraints on
deviations from general relativity due to an additional scalar degree of
freedom which are independent and competitive relative to constraints with
binary pulsar observations.




hector.okadadasilva@montana.edu
eXtreme Gravity Institute,
Department of Physics, Montana State University, Bozeman, Montana 59717 USA


nicolas.yunes@montana.edu
eXtreme Gravity Institute,
Department of Physics, Montana State University, Bozeman, Montana 59717 USA



Neutron star pulse profile observations as extreme gravity probes
    Nicolás Yunes
    December 30, 2023
=================================================================






§ INTRODUCTION

Extreme and mysterious, neutron stars are an
unavoidable consequence of general relativity at high enough
densities. With typical masses between 1 and 2 M_⊙,
but radii of only ∼ 11 km, their energy densities
can exceed nuclear saturation in their inner core. When this occurs,
matter can transmute into exotic forms that are impossible to
replicate in laboratories on Earth, and understanding the physics
of such extreme matter remains an open problem in
nuclear astrophysics.

Neutron stars are not just a laboratory for nuclear physics.
Their high densities also imply strong gravitational fields that exceed
those that can be probed in the Solar System by seven orders of
magnitude. Understanding the physics of such strong gravity objects
requires the use of a relativistic theory of gravity,
like general relativity <cit.>.

Neutron star observations can therefore provide invaluable
clues about both nuclear astrophysics and
relativistic gravity, as was achieved with the recent observation
of gravitational waves from the merger of two neutron stars <cit.>.

We here report on the first results of a new program to systematically
study how strong-field gravity can be probed with a particular set of neutron
star observations: pulse profiles emitted by radiating hot spots
on the star's surface <cit.>.

Neutron stars can have hot spots on their surface, i.e. regions
where the surface temperature is much higher than the average, due to
the impact of accreted material that is pulled in from a nearby, less dense star
or by localized heating due to magnetospheric currents (see
e.g. <cit.>).

These hot spots rotate rapidly with the star, emitting x-rays that
trace the spacetime geometry as they leave the star, producing a pulse profile
upon detection.

NASA's Neutron star Interior Composition Explorer (NICER)
is currently detecting these pulse profiles with unprecedented
(timing) resolution through rotation-resolved spectroscopy of thermal and
non-thermal emission <cit.>.

Can NICER constrain (or reveal) deviations from general relativity
in the strong-field regime of neutron stars?

The constraining power of such observations depend
sensitively on our ability to accurately model the pulse profile within and
outside general relativity. In <cit.> we presented a
complete toolkit to model x-ray pulse profiles in an entire class of
well-motivated modifications to general relativity: scalar-tensor
gravity <cit.>.

In this class, gravity is still described by the curvature of spacetime,
characterized by the spacetime metric as in general relativity,
but the metric is influenced by a scalar field that can be excited
in sufficiently dense environments, as in neutron stars.

Since the scalar field is usually negligible in the Solar System, scalar-tensor
theories have survived a plethora of experimental tests <cit.>,
while remaining prime candidates for tests with
neutron stars <cit.>.

The toolkit is fast, computationally efficient, covers a generic
family of scalar-tensor theories of gravity and it is ready to use for
probing strong gravity with NICER observations. The resulting pulse profiles
include Doppler shifts, relativistic aberration and time-delay effects,
thus extending the work of <cit.> to scalar-tensor theories
and the work of <cit.> to the required level of astrophysical
realism.

These effects are crucial for constructing sufficiently accurate
pulse-profile models, thus enabling, for the first time, a serious data
analysis investigation of the strength with which NICER observations
can probe strong gravity.

As a first step toward such an investigation,
we present an approximate, restricted Bayesian
calculation and find that a NICER observation consistent with general
relativity could allow for constraints on scalar tensor theories that can be
comparable to the most stringent constraints from binary pulsar
observations to date <cit.>.




§ LIGHT CURVE MODELING IN SCALAR-TENSOR GRAVITY


The x-ray photons emitted by hot spots at the surface of a neutron star trace
its exterior spacetime geometry. If the star is isolated, or if its companion
is far away and we focus on the neutron star's vicinity, the exterior spacetime
can be well-approximated by the Schwarzschild spacetime, provided the star
rotates much slower than the Kepler limit; this is indeed the case for some
targets of NICER, which have spin frequencies below 300 Hz.

In scalar-tensor gravity, the gravitational interaction is mediated by a
scalar field φ in addition to the usual metric tensor g_μν of
general relativity.

Here we consider a wide class of theories in which a massless scalar field is
minimally coupled to gravity, yet matter fields are coupled to the product of the
function A(φ) and the metric tensor g_μν. The function A determines the scalar-tensor model
under consideration.

The presence of the scalar field modifies the exterior spacetime of a static,
spherically symmetric star which is no longer the usual Schwarzschild spacetime
of general relativity, but instead is the Just
spacetime <cit.>. This spacetime depends not only on the
baryonic mass of the star m, but also on the “charge” or field strength q of the
scalar field φ.

The modifications relative to the Schwarzschild spacetime are
controlled by the scalar-charge-to-mass ratio
Q ≡ q c^2 / (G m) <cit.>, which is zero
in general relativity.

See <cit.> for a detailed discussion.


As a concrete example of scalar-tensor gravity we consider the two-parameter
theory of <cit.> in which

ln A(φ) ≡α_0 φ + β_0 φ^2 / 2.

The parameter α_0 is related to the Brans-Dicke
parameter ω_ BD via α_0^2 = 1/(2 ω_ BD + 3),
and it is stringently constrained to α_0 ≲ 3.4 × 10^-3 by
Shapiro time-delay measurements of the Cassini spacecraft <cit.>.

On the other hand, β_0 is unconstrained by Solar System experiments.

When β_0 ≳ -4.35, neutron stars have a small
Q ∝α_0 <cit.> and they are similar to those
of general relativity.

When β_0 ≲ -4.35 <cit.>,
regardless of the equation of state (EoS) <cit.>,
the theory admits a new branch of solutions, where stars with
Q ≈ O(1) are energetically favored over
their general relativity counterparts, due to a non-perturbative
effect called spontaneous scalarization[When the cosmological evolution of
the scalar field is taken into account, β < 0 violates current
Solar System constraints unless one fine-tunes the initial conditions
for the cosmological evolution. When β_0 > 0, these constraints
are satisfied <cit.> at the cost of disallowing
scalarization, unless β_0 ≳ 100 <cit.>.

We here study this theory with β < 0 as a toy-model for non-perturbative
strong-gravity effects, generically predicted in gravitational
theories with extra degrees of freedom (see e.g. <cit.>),
and because it allows for comparison with binary pulsar constraints <cit.>.
] <cit.>.

Numerically integrating the stellar structure equations in this theory,
we obtain m, q and A(φ), which completely determine the spacetime.

Figure <ref> shows that the scalar field
can significantly modify the stellar structure and thus the
spacetime around it.



With the spacetime geometry at hand, we can construct a
pulse profile model. In general relativity, the models used for
NICER data analysis are semi-analytic, as full null geodesic
ray tracing for the construction of millions of profiles in a Markov-Chain
Monte-Carlo exploration of the likelihood is computationally
prohibitive. A simple yet accurate model commonly used in general relativity is
Schwarzschild+Doppler <cit.>,
whose generalization to scalar-tensor theory leads to the “Just+Doppler"
model <cit.>, which includes Schwarzschild+Doppler
in the Q=0 limit.

As in the Schwarzschild+Doppler case, the Just+Doppler model
is constructed under some simplifications: the star is assumed spherical
(i.e. deformations due to rotation are ignored) and the effects of frame-dragging
in the photon's geodesic are neglected.

As in general relativity, these effects in scalar-tensor theory are
likely to have a negligible influence on the resulting pulse profile, as long as
the star rotates slowly (≲ 300 Hz).

This is because for slowly-rotating stars
we do not expect much of a deformation away from sphericity due to rotation.

Moreover, the frame-dragging of inertial frames in scalar-tensor gravity
is very similar to that of general relativity, justifying extrapolating
the conclusions from one theory to the other <cit.>.

For our analysis, we further assume that the hot spot is small (relative to the
size of the star) and isotropically radiating.

In particular, we assume that the spot has a small angular radius Δθ
and that it irradiates with a blackbody spectrum of constant k_ B T
measured in the frame comoving with the hot spot (k_ B being
Boltzmann's constant), while the rest of the star is dark.
As often done in the literature, we consider a single hot spot (see e.g. <cit.>)
and (for simplicity) neglect any other radiation coming from the background.
Finally, we assume that the observer collects photons in the
soft x-ray band in a single energy channel of 1 keV.

Although these simplifications are too rough to be implemented directly
in an actual analysis of NICER data <cit.>,
they are sufficient for a first data analysis
study, as shown through numerical simulations <cit.>
and analytical estimates <cit.>.

In scalar-tensor gravity, the model for the pulse profile depends not only on the
parameters

θ = {m, θ_ s, ι_ o, Δθ, D, k_ BT, f, EoS}

that describe the source, its geometry and the observer in general relativity, but
also on the theory parameters α_0 (we use log_10|α_0| in practice)
and β_0.

In this parameter space, m is the mass of the star, which given a specified
EoS determines its radius R, D is the distance to the source,
f is the star's rotation frequency,
θ_ s is the colatitude of the center of the hotpot and
ι_ o is the inclination of the line of sight, both angles
measured relative to the star's rotation axis.

Since we must first calculate neutron star models in scalar-tensor gravity,
instead of using m we use the central energy density
ε_15≡ε_ c× 10^-15 g/cm^-3
as our model parameter.

In summary, our full set of model parameters are
θ = {ε_15, θ_ s, ι_ o, Δθ, D, k_ BT, f, EoS, log_10|α_0|, β_0 }.




§ PROJECTED CONSTRAINTS FROM OBSERVATIONS



Given the model described above, let us now roughly estimate the accuracy to which
scalar-tensor theories could be constrained given a NICER observation
that is consistent with general relativity.

That is, we generate a synthetic signal (or injection) d
using the Schwarzschild+Doppler model in general relativity, and we
then attempt to extract it and estimate its parameters with the
Just+Doppler model.

A detailed analysis of the full likelihood is important, but beyond the scope of this
work if one is concerned only with forecasting preliminary constraints. If the forecasted
constraints were to reveal that NICER cannot place interesting bounds on scalar-tensor
gravity, then a more detailed analysis would not be necessary. However, in the opposite case,
one would conclude that NICER may be used as a laboratory to
test general relativity, and a more detailed data analysis study would be justified.
As we will show next, the latter is the case here.

The parameter space of the model is ten-dimensional and thus,
to arrive at a rough estimate on projected constraints, we
make a few simplifications that reduce this dimensionality.

We fix {θ_ s, ι_ o, Δθ, D, k_ BT, f, EoS}
at their injected values (indicated by asterisks). For the first five, we set

{θ_ s^∗, ι_ o^∗, Δθ^∗, D^∗,
k_ BT^∗} = {90^∘, 90^∘, 0.01^∘, 600 pc, 0.35 keV}.

Our choices for θ_ s^∗ and ι_ o^∗ are motivated by the
results of <cit.>, which showed that this is the best-case orientation for
determining the mass and radius of the source, with a 5% accuracy when the distance
to the star is also known (as assumed here).

The value of Δθ^∗ is chosen to enforce the point-like hot spot approximation,
while we chose typical illustrative values for D^∗ and k_ B T^∗.

Next, for EoS^∗
we consider both a soft (SLy4) and a stiff (APR) EoS in order to
test the variability of our constraints with the EoS. With an EoS^∗ chosen,
we choose ε^∗_15 to yield a neutron star with 1.93 M_⊙
as, e.g. is the case of the neutron star in the PSR J1614–2230 binary <cit.>.

This choice is motivated by the observation that high-mass stars are more susceptible
to the presence of the scalar field than low-mass ones (cf. Fig. <ref>).

We also consider two values of f^∗, 200 and 600 Hz, the latter
to magnify the contributions of Doppler and relativistic aberration effects on the
pulse profile[For these large rotation frequencies one should
include the influence of stellar oblateness in the model <cit.> using
e.g. the Oblate+Schwarzschild approximation <cit.>. As
shown in these papers, when θ_ s, ι_ o are close to the
equator, the effects of oblateness on the pulse profile are suppressed,and we expect
the same to be true in scalar-tensor gravity <cit.>.] <cit.>.

Finally, to generate our general relativistic injection, we
set log_10 |α^∗_0| = -5 and
β_0^∗ = 0[As noted earlier, these stars have Q ≈ 10^-5 and
are thus indistinguishable from neutron stars in general relativity as far as
NICER is concerned.].

With these simplifications, our parameter space becomes three-dimensional,
spanning only the central energy density ε_15 and the
theory parameters log_10|α_0| and β_0.

For ε_15, we assume uniform priors in the ranges
ε_15∈ [0.7, 3.0] (if EoS^∗ = SLy4) and
ε_15∈ [0.5, 2.4] (if EoS^∗ = APR), which
cover neutron star models with (roughly) minimum masses of 1 M  and
beyond the maximum mass for the corresponding EoS.

For β_0, we also use an uniform prior
β_0 ∈ [-7, 0], which includes the regime where
spontaneous scalarization happens.

The lower-bound already violates the best bounds from binary-pulsar
observations by more than 1σ <cit.>,
whereas the upper bound is chosen to include the limit of general relativity.

Finally, for log_10|α_0| we choose an uniform prior
log_10|α_0| ∈ [-5, -2]. The lower bound causes negligible
changes to the structure of neutron stars, whereas the upper bound
already violates the Cassini constraint by approximately one order of magnitude.



Having reduced our parameter space and chosen our priors, we calculate
the best-fit parameters by minimizing the relative
chi-squared between the injection and the model pulse profiles, scanning the
parameter space using 16 phase stamps over the course of one stellar revolution.
The standard deviation of the distribution σ_ε_15
is modeled as in <cit.> and chosen
to capture the optimistic 5% accuracy at which NICER can infer m and R.
From the chi-squared, we calculate the likelihood function and from it
we obtain the marginalized posterior distributions
p(ε_15|d), p(log_10|α_0| | d) and p(β_0|d).

Figure <ref> summarizes our results
for f^∗ = 600 Hz and EoS^∗ = SLy4.
The posterior of log_10|α_0| is essentially
flat, so little information is gained relative to the uniform prior.

The posterior of β_0 is
more interesting, clearly exhibiting a sharp decay that starts near
the scalarization threshold at β_0 ≈ -4.35.

When β_0 ≳ -4.35,
neutron stars are nearly identical to those of general relativity resulting in a
flat marginal posterior in this region.

However, when β_0 ≲ -4.35, the conditions for
spontaneous scalarization are realized, and, as
neutron stars deviate more and more from their general relativity
cousins as β_0 → -7, they result in poor models to recover the injection,
yielding an almost zero marginalized posterior when β_0 = -7.

The posterior for ε_15 peaks within ≲ 1% from the injected value ε^∗_15,
with a slow decay as ε_15→ 3. This is because
stars with the SLy4 EoS have a small variation in their radii and masses
past the maximum mass (≈ 2 M), which is close to the injected value
(1.93 M).

The conclusions above are robust to changes in the injected rotation frequency
and the EoS. When we set f^∗ to 200 Hz and the EoS^∗
to APR, the likelihood surface does not change considerably, although
the constraints on β are more stringent than when injecting SLy4,
as summarized in Fig. <ref>. The latter shows the
68% credible regions in the (|α_0|,β_0)-plane for both EoS injections,
which ought to be compared with the (1σ) constraints from Cassini and from
the absence of scalar dipole radiation from the PSR J1738-0333<cit.>,
PSR J0348-0432 <cit.> binaries <cit.>
and the hierarchical triple system PSR J0348-1715 <cit.>.

The contours reveal that for stiff equations of state (such as APR) the
constraints on scalar-tensor gravity (in the |α_0| → 0
portion of the parameter space) from NICER observations can, in principle,
be as strong as those obtained by binary pulsar observations[
A note of caution is necessary when comparing the constraints obtained here
and Refs. <cit.> shown in Fig. <ref>.
The constraints on β_0 depend critically on the choice of
the prior range and a naïve comparison between the three works is strictly
not correct. Moreover, the curves shown from <cit.> refer
specifically to EoS SLy4, whereas the EoS used in <cit.> (labeled “0.2” in <cit.>)
is considerably stiff, with a 2.6 M_⊙ maximum mass in general relativity.
].

Let us further compare and contrast the qualitative behavior of the constraint contours of Fig. <ref>.

As mentioned previously, β_0 is unconstrained from Solar System experiments, and
therefore, the Cassini contour is constant at a fixed |α_0|.

For binary pulsars, when β_0 > -4.35 neutron stars have small Q≈|α_0|,
and the modification to the binary's orbital time decay Ṗ_b scales with |α_0|^2.
But this modification is enhanced by two powers of the orbital velocity, which is of O(10^-2-10^-3),
because the scalar field generates dipole radiation, instead of the quadrupole radiation of general relativity.

The combination of these effects results in a modification to Ṗ_b that can be constrained
even for values of (α_0,β_0) beyond the threshold of scalarization and below the
Cassini bound.

Finally, in the NICER case, the modifications to the spacetime are controlled by Q, which is
small and roughly independent of α_0 inside our prior range and outside the threshold of scalarization
(i.e. outside the shaded region in Fig. <ref>). When scalarization does occur, however, Q acquires
a non-perturbative enhancement, that leads to a large deviation from general relativity, and thus to nearly vertical
constraint contours.






§ CONCLUSIONS

NICER allows for the observation of x-ray pulse profiles from rotating
neutron stars with unparalleled timing resolution. As the signal
travels from the star to the detector, it probes the spacetime in the
star's vicinity, opening the possibility of glimpsing (or
constraining) deviations from general relativity in the strong-gravity
regime.

Using the tools developed in <cit.>, we have
explored whether such observations can place constraints on a concrete
scalar-tensor gravity model which admits large deviations from
general relativity only in strong-gravity environments <cit.>.

As a proof-of-principle, we carried out a restricted likelihood analysis
that suggests that (in principle) observations carried by NICER can constrain the
parameter space of scalar-tensor gravity, severely restricting the domain in
which spontaneous scalarization occurs.

These projections are clearly preliminary but they do demonstrate,
for the first time, the potential of including pulse
profile observations as a new tool in the relativists' arsenal to perform strong-gravity tests.
This is of interest not just for the ongoing NICER mission, but also for planned,
future missions, such as the Enhanced X-ray Timing and Polarimetry
mission <cit.>.

Given these preliminary projections, there is now strong justification
to carry out a more detailed data analysis study that includes
a more careful modeling of the instruments's response, as well as the
covariances between all parameters in the model
to determine how much degeneracies impact our ability to test
general relativity.

The encouraging results found here motivate the consideration of other questions.
For instance, one could use both a general relativistic
and a scalar-tensor model to analyze real NICER data. Using Bayesian analysis,
one could then perform a model selection study to determine which theory is better
supported by the data through the Bayes factor.

Another important question is that of fundamental bias <cit.> in the
parameter estimation of the properties of neutron stars.
Pulse profile models assume a priori that general relativity is correct,
which may bias parameter estimation, if the true underlying
theory is not Einstein's.

As a consequence, the mass and radius of the source could be systematically biased with immediate
implications for EoS inference.

This question could be investigated through a Bayesian analysis that uses
a scalar-tensor theory pulse profile injection, but general relativity to model the data.


Still in the realm of synthetic data, a number of improvements can also be made.

First, we could include beaming of the emitted radiation
caused by Thomson scattering in the stellar atmosphere <cit.>.

Second, as radiation travels towards the observer, it interacts
with the star's accretion disk <cit.> and the
interstellar medium <cit.>.

Third, one could include background
radiation <cit.> and
generate synthetic data that is as realistic as possible.

This data would then have to be convolved with the detector's response matrix,
and combined with a Poisson sampling of the modeled pulse
profile to determine the chi-squared through a Poisson likelihood <cit.>.

With these ingredients taken into account, it would be interesting
to verify whether the constraints on scalar-tensor gravity obtained here are
robust to the assumptions we have made throughout the text.

Regardless of the simplifications made here, our work nonetheless puts forward the enticing
potential of probing strong-gravity with pulse profile observations.



§ DATA ANALYSIS


In this appendix we present a brief overview of the methods used to
obtain the projected constraints obtained in Sec. <ref>.



We call the signal measured during an observation the
synthetic injected signal, or (for brevity)
the injection F_ inj (θ^∗).

The pulse profile that we use to extract and characterize this observed pulse profile
is referred to as the model F_ mod(θ).

Both pulse profiles can be calculated following <cit.> once all
parameters
θ = {ε_15, θ_ s, ι_ o, Δθ, D, k_ BT, f, EoS, log_10|α_0|, β_0 }
have been specified.

As discussed in the main text, we consider the reduced model parameter space
obtained by fixing
θ_ fix = {θ_ s, ι_ o, Δθ, D, k_ BT, f, EoS}
to the injected values, leaving as variable model parameters
θ_ var = {ε_15, log_10 |α_0|, β_0 }

We calculate the best-fit parameter values by minimizing the
reduced chi-squared χ^2_ red between the injection and the model pulse profiles,
sampling over the model's variable parameters θ_ var.

The reduced chi-squared is defined as


    χ_ red^2 ≡
    
    
    1/N∑_i = 1^N[
    F_ mod(ϕ_i, θ_ fix, θ_ var) - F_ inj(ϕ_i, θ_ fix, θ^∗_ var)/σ(ϕ_i)]^2 ,


where the summation is over N time-stamps during the course of one revolution
of the star. We normalize the phase over a revolution such that ϕ_i ∈ [0, 1],
and use N = 16 time stamps.

The standard deviation of the distribution σ is modeled as


    σ(ϕ_i) =
    σ_ε_15(ϕ_i) +
    σ_log_10|α_0|(ϕ_i) +
    σ_β_0(ϕ_i) ,


where σ_ε_15, σ_log_10|α_0| and σ_β_0 are the
standard deviations on each injected parameter
θ^∗_ var = {ε_15^∗, log_10 |α_0^∗|, β_0^∗}.

Since our injection is assumed to be consistent with general relativity, we set
σ_log_10|α_0| = σ_β_0 = 0,

while the standard deviation σ_ε_15 is calculated by <cit.>


    σ_ε_15(ϕ_i)    = 1/2|
    F_ inj(ϕ_i, θ_ fix, {ε_15^∗ + δε_15^+, log_10 |α_0^∗|, β_0^∗})  .
    
         .
        -    F_ inj(ϕ_i, θ_ fix, {ε_15^∗ - δε_15^-, log_10 |α_0^∗|, β_0^∗})
    | ,


where (as explained in the main text)
log_10 |α_0^∗| = -5 and β_0^∗ = 0,
while δε_15^± remain to be specified.

To obtain them, we start by assuming that m^∗ = 1.93 M_⊙ as
discussed in the main text.

Then, given an EoS, we calculate the
central energy density ε_15^∗ for which this mass is
obtained and record the corresponding radius R^∗.

Next, we draw
a box in the mass-radius plane (in general relativity) centered at (m^∗, R^∗)
with width spanning R^∗(1 ± 0.05) and height spanning m^∗(1 ± 0.05).

The mass-radius curve intersects this box twice and these two points then determine
δε_15^±.

Once the reduced chi-squared is obtained, we calculate the likelihood


    L(θ_ var) = exp(- χ_ red^2 / 2) .


When this is evaluated over the whole prior domain, we obtain the likelihood distribution
shown in the off-diagonal panels of Fig. <ref>. Upon marginalization of the likelihood,
we obtain the diagonal panels of the same figure.

To obtain the constraints shown in Fig. <ref>, we first sort
our likelihood data in decreasing order of values of L(θ_ var).

Then, we sum all the values of
L (and record the corresponding θ_ var) until 68% of the
sum over all L-values is reached.

The values of log_10|α_0| and β_0 obtained from this procedure
result in the dotted and solid lines in Fig. <ref>.






§ ACKNOWLEDGMENTS


We thank Anne Archibald, Justin Alsing, Alejandro Cárdenas-Avendaño,
Neil Cornish, Cole Miller, Sharon Morsink, George Pappas and
Hajime Sotani for valuable discussions.

We also thank David Anderson for calculating and sharing with us the
binary pulsar constraints.

Finally, we thank the anonymous referees for pertinent comments on
our work.

This work was supported by NASA grants NNX16AB98G and 80NSSC17M0041.

Computational efforts were performed on the Hyalite High Performance
Computing System, operated and supported by  University Information Technology
Research Cyberinfrastructure at Montana State University.

apsrev4-1


