Warp Factory: A Numerical Toolkit for the Analysis and Optimization of Warp Drive Geometries

Abstract

The last few decades of warp drive research have focused on analytic methods to explore warp solutions to Einstein's field equations. These analytic solutions tend to favor simple metric forms which are easier to analyze but limit the space of exploration. In addition, all solutions to date have involved unphysical qualities, such as negative energy, violation of energy conditions, or enormous energy requirements. In an effort to explore the space of physically meaningful warp drives, the Advanced Propulsion Laboratory (APL) at Applied Physics has developed Warp Factory, a toolkit written in MATLAB for numerically analyzing and optimizing warp drive geometries. Warp Factory consists of a series of three primary modules: the solver, the analyzer, and the optimizer. Together, these modules allow users to solve the Einstein field equations, compute energy conditions and scalars, and perturbatively optimize general metrics. Finally, the toolkit offers insightful 2D and 3D visualizations of general metrics and stress-energy tensors. The methods used in Warp Factory, along with their application in evaluating and optimizing common metrics, are discussed. With Warp Factory, APL hopes to accelerate warp research and bring us one step closer to physical and realizable warp drives.

I. Nomenclature

• α — lapse rate
• βⁱ — shift vector
• γᵢ — spatial terms
• κᵢ — shear terms
• Γ^α_μν — Christoffel symbol
• ρ_E — energy density
• ρ_NEC — null energy condition
• ν_WEC — weak energy condition
• k^a — null vector field
• X^a — timelike vector field
• X — spacetime 4-vector
• p_σ — momentum density
• S^αβ — pressure / shear tensor elements
• η_μν — Minkowski metric
• M^μ_ν — frame transformation matrix
• g_μν — metric tensor
• T^μν — stress-energy tensor
• R_μν — Ricci tensor
• R — Ricci scalar

II. Introduction

Since the publication of Introducing Physical Warp Drives the Advanced Propulsion Laboratory (APL) at Applied Physics has continued its exploration of physically meaningful warp drives. To date, the efforts in warp research have been constrained to an analytical approach, which has often provided a limited, and in some cases, inaccurate, assessment of the viability of proposed metrics. To address the challenges of analyzing warp drives, the APL team has developed a numerical analysis toolkit called Warp Factory, which provides a fully generalized solver of the Einstein Field Equations to evaluate the stress-energy tensor along with important scalars and energy conditions. In addition, Warp Factory also provides a novel optimization approach that uses a perturbation-based technique to generate new metrics that can optimize the physicality conditions of a starting warp metric.

In the following paper, an overview of the Warp Factory codebase and its analysis of common warp metrics will be presented. In Section III, the concepts and general equations of General Relativity (GR) are presented. In Section IV, the modules of the code are described along with its implementation. In Section V, popular and recently proposed metrics are numerically evaluated and analyzed. In Section VI, the optimization algorithms are demonstrated. Finally, in Section VII, the progress and development of the Warp Factory codebase along with the future work planned for Warp Factory and the APL team will be summarized.

III. Background

A. Metric and Stress-Energy Tensor

GR is a geometric theory that describes the effects of gravity as a curvature of spacetime. It provides the mathematical machinery to connect the curvature of spacetime with the energy, pressures, and stresses needed to create that curvature. At the center of this theory are the Einstein Field Equations (EFE), which connect the metric tensor g_μν to the stress-energy tensor T_μν via the Ricci curvature R_μν and Ricci scalar R:

R_μν − ½ R g_μν = (8πG/c⁴) T_μν (1)

where the Ricci curvature R_μν is built, in Eq. (2), from the derivatives of the Christoffel symbols with respect to the coordinates, minus the corresponding derivative term with the indices contracted, plus the quadratic combinations of pairs of Christoffel symbols, summed over the repeated indices — which is constructed using the Christoffel symbols Γ^γ_αβ:

Γ^γ_αβ = ½ Σ_σ ( ∂g_βσ/∂x^α + ∂g_ασ/∂x^β − ∂g_αβ/∂x^σ ) g^γσ (3)

The Ricci Scalar R is:

R = Σ_μ Σ_ν g^μν R_μν (4)

Henceforth, all Greek indices will range from 0 to 3, whereas Latin indices range from 1 to 3. Summations will be implied following Einstein notation.

The metric tensor is a 4x4 matrix whose time-time entry is −α² plus the contraction of the shift vector with itself, whose time-space and space-time entries are the shift vector components, and whose space-space block is the spatial metric γ_ij (5)

Here the metric terms are expressed following the 3+1 formalism and using the minus-plus-plus-plus metric signature. It should be noted that the metric is symmetric, hence the duplicated indices across the diagonal. In 3+1 the α is the lapse rate which describes the passage of time, βⁱ is the shift vector describing spatial motion, and γ_ij are the spatial terms which describe spatial shearing when the indices differ and the expansion or contraction of space when they are equal. Likewise, the stress-energy tensor is also a 4x4 matrix, with a time-time entry T00, time-space entries T0i and T0j, and a space-space block Tij (6)

Like the metric tensor, this is also symmetric. The T00 represents the energy density, the T0i is the momentum density of the associated energy, Tii are the pressures, and Tij are the shear stresses, as observed by "at rest" observers.

When constructing the metric, both the coordinate system used to express the metric and the physical frame are important to consider. They are defined as follows:

• Coordinate System: The spacetime grid selected to represent the space; examples are Cartesian and spherical coordinate systems.
• Frame: In addition to the coordinate system, the frame of observation defines an additional meaning to the coordinates; an example is a comoving frame where time-varying elements are rolled into the definition of the coordinates themselves. Frames define what observers see from their perspective.

For warp metrics, the selection of the coordinate system can reduce the amount of computation by exploiting relevant symmetries. Since warp drives have a direction of travel it is natural to conduct analysis in an axisymmetric coordinate system, which exploits a φ symmetry in the metric around the direction of travel. Additionally, performing analysis in a comoving frame can reduce the need to consider time derivatives for constant velocity warp drives when solving the EFE.

When interpreting the metric and stress-energy tensors it is important to consider what frame these are defined within, as this will impact the meaning of the terms. Different observers in general will see different values. Frame transformations can be done directly on the stress-energy and metric tensor to transform between observers:

g^μν = M^μ_α M^ν_β g^αβ, and T^μν = M^μ_α M^ν_β T^αβ (7)

where M^μ_ν is a transformation matrix. This is defined by a required transformation to shift the metric from one frame into another, which is then also applied to the stress-energy tensor as shown in Equation 7.

Finally, all observable elements are found by taking the stress-energy tensor and contracting it along a defined observer. The contraction process determines how a given measurable quantity is determined from the tensor as measured by the observer. For any observer, the measurable energy density might involve combinations of different tensor terms when compared from one frame to another. The tensor is not invariant between frames. In general, the contraction for energy density ρ, momentum density pⁱ and stresses S_ij are given by:

ρ_E = T_μν n^μ n^ν (8)

p_σ = − T_μν n^μ γ^ν_σ (9)

S_ij = T_μν γ^μ_i γ^ν_j (10)

where n^μ and γ^μ_ν are the projection vectors of the observer. Other scalar values will be discussed in Section IV.B.4.

B. Warp Structure

In this analysis of warp drives, some common terminology will be used in describing their geometry and features throughout this paper. These important terms are described below:

• Passenger Volume: The inner volume where the "passengers" or cargo of the warp drive exists, shifted along by the warp bubble around it.
• Warp Bubble / Bubble Region: The non-flat region which transitions the passenger volume to the world boundary and where the required matter-energy of the drive exists.
• Observer: Someone or something moving along a timelike future-directed worldline who can make local measurements to the world around them. These observers are expressed as having 4-vectors with their time and spatial terms restricted to values satisfying a timelike condition. A special observer, often used in the analysis of metrics and tensors, is the Eulerian observer. This is an observer that resides within a frame locally at rest. In spacetime, this means an observer not moving through space but progressing forward in time. (A timelike vector is defined as a 4-vector X whose inner product with itself, formed with the metric, is less than zero.)

The APL team's working definition of a warp drive is a metric structure that has and/or creates the following effects on observers in spacetime:

1. Drags passengers along with the moving warp bubble. This is a unique feature of a warp metric, as it must provide an effect on observers which is different from having a force applied directly.
2. Flat or constant spacetime within the passenger volume. The inner region should be flat according to any and all observers situated within the warp bubble.

There are a few other unique properties that warp drives can provide compared to typical modes of transportation. These are useful and desirable features but not requirements for a warp drive. The features of interest are described below:

• Contactless acceleration of passengers. Does the warp structure give an effective acceleration to observers without a force imparted to the passengers? This could allow very high accelerations to be achieved without stress on the passengers or cargo inside.
• No momentum flux to infinity. Does the drive shed its matter-energy to accelerate observers? Likely a change in the energy is required to accelerate the bubble, but whether this is dominated by gravitational waves or classical ejection of matter, like a rocket, is an open topic.

C. Physicality

A meaningful warp metric should be something which is physical and hopefully constructable. The physical nature of a warp drive can be expressed as a set of conditions applied to the stress-energy tensor. These conditions are relationships between what any observer would see in the energy density, momentum, pressures, and stresses that align with what would be physical. These energy conditions are the following:

• Non-negative energy in the Eulerian frame. Energy density should be positive in the locally flat frame of an Eulerian observer.
• Satisfies the Null Energy Condition (NEC) everywhere. The NEC expresses the observed mass-energy density for null (light-ray) observers. For this to be physical, mass-energy density must be non-negative as viewed by any null or lightlike observer.
• Satisfies the Weak Energy Condition (WEC) everywhere. The WEC expresses the observed mass-energy density for timelike (matter) observers. For this to be physical, mass-energy density must be non-negative as viewed by any timelike observer.
• Satisfies the Dominant Energy Condition (DEC) everywhere. The DEC expresses the velocity of observed matter flow. For this to be physical, the matter should not be seen flowing faster than light by any timelike observer.
• Satisfies the Strong Energy Condition (SEC) everywhere. The SEC expresses the tidal effect of gravity acting on observers. For this to be physical, the gravitational tidal effect should be non-negative, i.e., matter always gravitates towards other matter.

In addition to constraints on the matter-energy conditions, to have a realizable warp metric further conditions should be applied based on the practical limits to construct the warp drive:

• Known stress-energy tensor form. The required elements of the stress-energy tensor should follow those of physically known systems. An example of an unknown stress-energy tensor is one that has pressure terms but no associated energy density.
• Uses reasonable mass. Any warp bubble should require a practically attainable amount of matter-energy.

Finally, while not an immediate requirement for a physical warp drive, it is of interest to consider the properties when the velocity of the system goes superluminal, namely: are physicality and causality maintained during superluminal motion?

In this paper, the focus will be on evaluating a few of the energy conditions, limited to only the Eulerian energy, NEC, and WEC.

IV. Toolkit Structure

The Warp Factory code is built using MATLAB and consists of several functions, organized into modules, that perform different aspects of the overall analysis, visualization, and optimization. In Figure 1, a basic structure of the Warp Factory codebase is shown. There are three main modules within Warp Factory: the Solver, Analyzer, and Optimizer.

A. Solver Module

The solver module handles the solution methods to the EFE. The EFE are solved using a second and fourth-order central finite difference scheme, employed on the spacetime grid. The solver function can use either the CPU or GPU of the host machine, with the selection depending on the grid size for best performance. For a basic Cartesian grid, the finite difference solvers assume the derivatives at the bounds are zero. To extend this solver to other coordinate systems, padding points are required for the correct computation of derivatives with radial coordinates or periodic boundaries. At this time, Warp Factory has focused on implementing only the Cartesian and cylindrical coordinate systems. For a 2D grid in cylindrical or axisymmetric coordinate frames, the addition of "ghost points" are required for the computation of derivatives around r equal to zero. These ghost points mirror the real points across r so that the minimum grid points are provided for the finite difference method at r equal to zero. The mirroring process also depends on the metric component, with some terms being mirrored symmetrically, whereas other terms are mirrored anti-symmetrically. If the size of the spacetime is defined to be unity along any dimension (for example, the time or a symmetric direction like φ), the solver assumes the derivatives along those dimensions are zero.

Fig. 1 — Diagram of the major functions and structure of the toolkit. Note that some blocks do reference each other outside of this general hierarchy. Also not shown are smaller functions and applications for plotting data.

B. Analyzer Module

The analyzer module takes in the metric and stress-energy tensor and generates observables for the user. This module also conducts the relevant transformations and contractions to generate these values. Two kinds of transformations can be performed: frame or coordinate transformations.

1. Frame Transformations

Warp Factory currently implements a frame transformation solver to convert the stress-energy tensor into an Eulerian frame from any arbitrary input metric. The Eulerian frame is defined as the frame which has a flat Minkowski metric:

η_μν = diag(−1, x1, x2, x3) (11)

The coordinate system impacts the frame transformation. For Cartesian, x = [1, 1, 1], but for cylindrical x = [1, r², 1]. (Note that typically the Minkowski metric is the same in covariant and contravariant forms, but this is only the case in Cartesian coordinates and not the case for a general coordinate system.) The transformation matrix solves, at all points in spacetime, the matrix equation:

Mᵀ g M = η (12)

where M is a lower triangular matrix. This transformation has a general symbolic solution, albeit a complicated one, which is computed for all points. With M defined, the stress-energy tensor is transformed using Equation 7. The motivation for this transformation is so that further point-wise analysis on the energy tensor components can be conducted within a well-defined frame, such as when determining energy condition violations which require defining a vector field.

In addition to Eulerian transformation, shifting the metric into the comoving frame saves computation for constant velocity warp drives by removing the need for derivatives along t. The transformation acts on the shift vector:

β′ = β − v (13)

This results in shifting the frame of analysis to that of a frame comoving with the passenger volume and is typically centered in the world volume, with the world boundary moving at minus v. In the comoving frame, there are no coordinates of the warp drive which need to vary with time, thus all time-derivative terms are removed.

2. Coordinate Transformations

Transformations between the 2D and 3D coordinates can be performed directly on the metric components, mapping the metric tensor from one coordinate to another. This mapping solves the transformations of the line element using standard coordinate relations. When mapping from axisymmetric 2D to Cartesian 3D, points beyond the radius boundary of the original 2D system are set to a specified boundary condition, due to a lack of data to extrapolate. To improve the accuracy of the translation, a 3rd-order Lagrange polynomial interpolation is used. These transformations also are conducted on the vector fields as needed and in the definition of the Minkowski metric, whose diagonal elements change depending on the coordinate system.

Contravariant and covariant transformations are done using the metric tensor acting on vectors or tensors as an index-raising, lowering, or mixing operation:

U_ν = g_μν U^μ, and U^ν = g^μν U_μ (14)

T_μν = g_νβ g_μα T^αβ, and T^μν = g^νβ g^μα T_αβ (15)

T^μ_ν = g_να T^μα (16)

The metric tensor g_μν is transformed into covariant or contravariant forms using matrix inversion at each point.

3. Energy Conditions

The primary mathematical expression of physicality to a warp metric is in the satisfaction of energy conditions. There are four main conditions used in the literature: null, weak, dominant, and strong. Currently, Warp Factory only implements analysis on the null and weak, which are described below.

The Null Energy Condition (NEC) is expressed at a given point in spacetime as the requirement that ν_NEC(X), formed by contracting T_ab(X) twice with a null vector k, is greater than or equal to zero (17)

where T_ab is the energy tensor transformed to an orthonormal frame at each point in spacetime X and for a future-pointing null vector field k defined (in flat space) as:

k = (t, r) (18)

where the null vector field satisfies the condition that the inner product of k with itself is zero (19)

(In GR the inner product is generalized as the metric contracted with the two vectors, which differs from 3-vector inner products.)

The value of t is fixed to unity and all spatial vectors in the field satisfy the condition that the magnitude of r equals one for a specified number of vectors. The density of vectors will define the resolution of finding the maximum null violation amount. Each spatial vector is sampled by generating a grid of uniformly distributed points on a sphere of unity radius. This process is done for each spacetime point with a map of the maximum violation returned. While this is often considered as a binary check on the condition that ν_NEC(X) is non-negative across all spacetime, the magnitude of violation can be used as a fitness parameter itself during optimization.

The Weak Energy Condition (WEC) is calculated in a similar manner, but uses a timelike vector field: ρ_WEC(X), formed by contracting T_ab(X) twice with the timelike vector X, must be greater than or equal to zero (20)

where X is defined similarly to k but as a timelike vector field satisfying the condition that its inner product with itself is less than zero (21)

The spatial terms of r in the vector field X are generated evenly across the surface of a sphere as before, but also in a series of shells with decreasing magnitudes of r to 0. These shells are evaluated and the maximum violation is returned for all directions and all shells. It is also important to note that NEC violation implies WEC violation.

4. Scalars

The properties of the metric can be explored by comparing coordinate invariant quantities, specifically the expansion and shear scalars. Each of these scalars is built from the metric tensor and evaluated at each point in spacetime using projections of the metric on a timelike vector field U^μ, defined by the lapse and shift vectors as:

U^μ = (1/α)(1, −βⁱ), and U_ν = g_μν U^μ (22)

(This is just one example of a timelike vector field of interest; these quantities can be computed for any timelike vector field.)

The two scalars of interest are expansion and shear.

The expansion scalar θ is defined as:

θ = g^μν θ_μν (23)

where θ_μν is the stress tensor found from the projection tensor P_μα acting on the observer vector field U as:

θ_μν = P_μα P_ν^β ∇_(α U_β) (24)

and the projection tensor is given by:

P_μν = g_μν + U_μ U_ν (25)

The shear scalar σ² is defined as:

σ² = ½ σ_μν σ^μν (26)

where σ_μν is the shear tensor found from the stress tensor θ_μν, the projection tensor P_μν and the expansion scalar θ as:

σ_μν = θ_μν − (θ/3) P_μν (27)

These two scalars can be interpreted as invariant quantities signifying certain effects the metric tensor has on the spacetime for the observer U. Specifically, they define the following:

• Expansion Scalar: The expansion scalar can be thought of as a quantity that describes the change in volume of an object.
• Shear Scalar: The shear scalar can be thought of as a quantity that describes the relative expansion or contraction along any dimensions while preserving the volume of an object.

The most common scalar discussed in the literature is expansion. For the Alcubierre metric, this exhibits a contraction in front of the bubble and an expansion behind it. However, it should be noted that the presence of expansion, or any of these scalars being non-zero, is not fundamentally required to make a warp drive.

5. Momentum Flow

Most warp metrics generate terms beyond just the T00, including pressures, shears, and momentum. The momentum describes the motion of the energy associated with the warp metric. To understand the flow of matter in the bubble a map of flowlines is constructed. This is done by building a velocity vector field comprised of the momentum flux terms pⁱ across each spacetime point:

Ω(X) = ( p¹(X), p²(X), p³(X) ) (28)

where T^μν is first transformed into the Eulerian frame and contracted to build pⁱ. The flowline is generated using a "virtual" timestep from a starting grid of points, which are pushed along by a velocity imparted by the momentum vector field. Note that this is not a trajectory of matter moving along a geodesic, but rather a visualization that captures the snapshot of the momentum density structure in the warp bubble as generated by the metric tensor. The virtual time steps are not real-time evaluations of the dynamical spacetime across t. Momentum flowlines are currently only generated for metrics in Cartesian coordinates, with coordinate transformations performed as needed. Since this flowline is computed on a finite grid, a trilinear interpolation is used on the momentum vector field data for interpolation between points.

C. Optimization Module

The complexity of solving the EFE presents a limitation to any analytical approach in exploring the space of complicated warp metrics. One method to dynamically generate new metrics is through a simple Machine Learning (ML) approach. This technique will perform perturbations to the metric, compute their impact on the stress-energy tensor, and then calculate the resulting energy conditions as a fitness measure on improving the physicality of a warp drive. Each perturbation is selected based on its improvement to the fitness.

The ML approach is implemented as a Monte Carlo Gaussian Perturbation within a Genetic Algorithm. To allow for general perturbations that will not break the metric, a new formalism is developed called 6+4. Under 6+4, the usual 3+1 approach is extended with modifications in order to preserve the metric signature minus-plus-plus-plus under arbitrary perturbations to the underlying functions which generate it. The form of this metric, given in Eq. (29), is a 4x4 array whose leading diagonal carries g00, g11, g22 and g33, whose first row and column carry the shift components β1, β2 and β3, and whose remaining off-diagonal entries carry the shear terms κ1, κ2 and κ3, where each of the diagonal terms is built from the ten independent 6+4 functions βᵢ, κᵢ, γᵢ and α as:

g33 = γ3² (30)

g22 = γ2² + (κ3/γ3)² (31)

and, in Eqs. (32) and (33), g11 and g00 are built in the same nested way from γ1, γ2, γ3, the shear terms κ1, κ2 and κ3, the shift terms β1, β2 and β3, and the lapse α, each successive term subtracting the contributions already carried by the terms below it in the sequence.

During the optimization, changes are made to each of the 6+4 functions and the full metric is then constructed using the equations above. The modifications made to enable 6+4 to be robust in the face of perturbations do not share the same physical interpretations which drive the formalism of 3+1 for a generic metric. While this is important for the perturbation code, the use of 3+1 decomposition in interpreting results should still be used.

1. Monte Carlo Gaussian Perturbation Method

The core element of the ML approach is the Monte Carlo Gaussian Perturbation (MCGP) technique. In MCGP each perturbation is applied, centered on points selected in the spacetime, for each function in the 6+4 formalism. These perturbations are Gaussian functions that can vary in both amplitude and σ. A minimum σ size bound is used to avoid sharp changes, which can exploit issues with the numerical grid resolution in the finite difference method. This process is inherently time-consuming as the perturbations require computation of the EFE and evaluation of energy conditions to determine their impact on fitness. This is improved by doing perturbations to parts of the space and applying symmetry where possible. The selection of which points to perturb are also selectively chosen and currently use the NEC violation map to focus changes around regions that need the most adjustment. The method uses a probability of selection which scales the violation logarithmically: in Eq. (34), the selection probability at a point is an amplitude a multiplied by the exponential of the base-ten logarithm of the null-energy violation at that point.

The amplitude of the perturbation is expressed as a fraction of the local metric value to enforce perturbations that are scaled appropriately.

The usage of symmetry starts by running the MCGP only on 2D axisymmetric grids. In addition to using a symmetric coordinate, perturbations in the ρ, z space are also mirrored around the midpoint along the z-axis (direction of travel), with most terms mirrored symmetrically, and the β2, κ1, and κ3 mirrored anti-symmetrically. The MCGP cycles through making perturbations to each metric function. Which metric functions are allowed to vary are user-selectable.

The acceptance or rejection of perturbations is based on improvement of fitness. To decrease the chance of getting stuck at a local minimum, the acceptance of the perturbations allows for a certain decrease in the fitness to be accepted, within a user-selectable margin. This tends to allow better improvement in longer runs. The boundaries of the world and passenger volume are held fixed during the MCGP.

2. Genetic Algorithm

A key feature of the MCGP is the randomness in its perturbations. To improve both the effectiveness and speed of execution, a genetic algorithm is used. This algorithm runs the MCGP on several mutated metrics in parallel. A mutated metric is one that has random changes applied to it. From these parallel runs or populations, the best result can then be selected and further mutated again and run in parallel across several generations. The process flow for this is shown at a high level in Figure 2. The mutation process adds Gaussian perturbations to a starting metric regardless of whether it improves or worsens the fitness. The addition of starting mutations allows even more robustness in the face of local minima.

Fig. 2 — Selection and mutation flow using the genetic optimization algorithm.

V. Metric Evaluations

Using Warp Factory, common and new warp metrics discussed in the literature can be evaluated from a numerical perspective, showcasing them in greater detail.

A. Alcubierre Metric (AM)

The most famous and very first warp metric was proposed by Alcubierre in 1994. This metric uses a flat passenger volume with a constant shift vector which connects to the boundary by a hyperbolic tangent function. The AM is shifted into a comoving frame, which simply inverts the shift vector to be plus v at the world boundary and zero inside the passenger volume; in this frame time derivatives can be ignored.

Using Warp Factory, we construct the AM in 3D coordinates and solve the EFE. The shift vector is shown in Figure 3, the Eulerian energy density in Figure 4, the NEC and WEC in Figure 5, the full stress-energy components in Figure 6, and finally the scalars in Figure 7.

Fig. 3 — Metric functions for the AM. Only the shift vector along Z is modified from Minkowski. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 4 — Eulerian energy density for the AM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

As identified by Alcubierre himself, the AM metric violates both the NEC and WEC within the bubble region due to its requirement of negative energy density; this result is confirmed here using the code.

More generally, the AM metric falls into what is commonly referenced as a Natário class warp metric which is built only using shift vector components: in Eq. (35) the metric takes the form of a 4x4 array whose time-time entry is minus the quantity one minus v squared, whose time-space and space-time entries are minus the components of v, and whose space-space block is the Kronecker delta,

where v follows the passenger volume of choice.

With Warp Factory, different and more complicated passenger volumes can be explored to see how the topology of the inner volume impacts the energy conditions. For this analysis, the passenger volume shift velocity is set along Z and is smoothly connected to the boundary using a Gaussian smoothed falloff over a specified radius R. The Eulerian energy and NEC are compared between a pancake, cigar, cylinder, bicone, and teardrop passenger volume shapes to that of a standard spherical volume. All shapes maintain a constant total passenger volume for an equal comparison.

Fig. 5 — NEC and WEC for the AM. Direction of motion is along +Z. Plotted are the most violating values across all observers. Cross-section is plotted for y = 0.

Fig. 6 — The stress-energy tensor, in the Eulerian frame, for the AM. The T00 component is shown in Figure 4. Direction of motion is along +Z. Cross-section is plotted for y = 0.

The shift vectors for each shape are shown in Figure 8. The energy density and NEC are shown in Figures 9 and 10. For each passenger volume, the integrated Eulerian energy density and NEC violations are used to compare how well each transport the same passenger volume at the same velocity. The results of this are shown in Table 1. The change in passenger volume shape does see an improvement in the amount of energy required and a reduction in the NEC violation, with the pancake and teardrop shapes performing the best. However, these improvements are all under an order of magnitude, which means the selection of passenger volume is not the primary factor in either physicality nor massively reducing the needed energy. In Figure 11, the impact of the different passenger volume shapes on the density and integrated Eulerian energy and NEC violations are shown.

Fig. 7 — Expansion and shear scalars for the AM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 8 — Shift vector along the Z axis for various passenger volume shapes. They are defined in a comoving frame. Cross-section is plotted for y = 0.

Fig. 9 — Eulerian energy density for each of the passenger volumes. Cross-section is plotted for y = 0.

Fig. 10 — NEC for each of the passenger volumes. Cross-section is plotted for y = 0.

Fig. 11 — Scatter plots of the Eulerian energy density, Eulerian total energy, total NEC violation, and NEC violation density relationships for the various passenger volumes.

Table 1 — Evaluation results for the different passenger volumes. The total Eulerian energy is a pointwise sum of the locally transformed stress-energy tensor T00 component across each point in spacetime, making this value a representation of the sum total of each perceived energy, for each observer, at each point. This value is not the total integrated energy as viewed from a single observer across spacetime.

| Parameter | Spherical | Pancake | Cigar | Cylindrical | Bicone | Teardrop | |---|---|---|---|---|---|---| | Total Eulerian Energy [J] | −4.2 × 10⁴² | −1.9 × 10⁴² | −6.3 × 10⁴² | −3.1 × 10⁴² | −3.8 × 10⁴² | −2.1 × 10⁴² | | Max Eulerian Energy Density [J/m³] | −2.6 × 10³⁸ | −2.0 × 10³⁸ | −2.6 × 10³⁸ | −2.6 × 10³⁸ | −2.0 × 10³⁸ | −0.9 × 10³⁸ | | Total NEC Violation | −3.0 × 10⁴⁴ | −2.0 × 10⁴⁴ | −4.0 × 10⁴⁴ | −2.5 × 10⁴⁴ | −3.0 × 10⁴⁴ | −2.4 × 10⁴⁴ | | Max NEC Violation Density | −1.0 × 10⁴⁰ | −0.7 × 10⁴⁰ | −1.2 × 10⁴⁰ | −0.9 × 10⁴⁰ | −3.0 × 10⁴⁰ | −2.3 × 10⁴⁰ |

B. Van Den Broeck Metric (VDM)

Following Alcubierre's publication, Van Den Broeck built a modification to the AM in 1999 by using concentric regions which varied the shift and expansion components in an effort to reduce the negative energy required via the expansion of the passenger volume. The metric components cross-section is shown in Figure 12 demonstrating these regions in the Z shift vector value. The Eulerian energy density has a similar shape to that of Alcubierre but has a positive region on the inside through the transition into the passenger volume with its expansion. This is shown in Figure 13. Not surprisingly, this metric has WEC violations similar to Alcubierre, shown in Figure 14. The existence of negative energy in the Eulerian frame is sufficient to know that a violation exists but the reason for different rings of violation is better understood when considering the stress-energy tensor, shown in Figure 15. The interaction of the pressure and momentum density terms when the energy density is small causes the main violation regions.

Fig. 12 — Metric functions for the VDM. The VDM uses both the shift vector along Z and passenger volume expansion terms along X, Y, and Z. Cross-section is plotted for y = 0.

Finally, the expansion and shear scalars are shown in Figure 16. Each of the transition regions between the passenger volume and bubble has expansion and shear present, in a form similar to the AM.

Ultimately, the VDM does not provide an approach that resolves WEC violation, but it does provide a novel method to make a warp drive more efficient, as the expansion terms allow more passenger volume per required energy.

Fig. 13 — Eulerian energy density for the VDM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 14 — NEC and WEC for the VDM. Direction of motion is along +Z. Plotted are the most violating values across all observers. When only positive energy density is seen across all null or timelike observers, then the minimum positive value is shown. Cross-section is plotted for y = 0.

Fig. 15 — Eulerian stress-energy tensor for the VDM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 16 — The expansion and shear scalars for the VDM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

C. Bobrick-Martire Modified Time Metric (MTM)

In Introducing Physical Warp Drives, Bobrick and Martire consider several modifications to the AM, one of these being the addition of a changing lapse rate (see section 4.5 of that paper). The metric functions for this kind of drive are shown in Figure 17. The impact of a changing lapse is a reduction in the passage of time for passengers inside the drive as compared to an outside observer. The Eulerian energy density resulting from an addition of changing passenger lapse is remarkably similar to Alcubierre, with the energy density shown in Figure 18. However, when considering the energy conditions, a much different result is found, shown in Figure 19. For both NEC and WEC there is no longer any violation inside the inner region, with the NEC only seeing positive values within that region. Looking at the stress-energy tensor in Figure 20, positive pressures are seen on the inner region of the bubble, as opposed to the negative values seen in the AM. Finally, the expansion and shear scalars are shown in Figure 21. These share the same shape and values as seen in the AM.

While the addition of a changing lapse to the metric is not a fix for all violations of the WEC, it does offer an improvement over the AM by removing violations in the inner region of the bubble.

Fig. 17 — Metric functions for the MTM. Only the lapse and shift vector along Z have non-Minkowski values. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 18 — Eulerian energy density for the MTM. Cross-section is plotted for y = 0.

Fig. 19 — NEC and WEC for the MTM. Direction of motion is along +Z. Plotted are the minimum values across all observers. When only positive energy density is seen across all null or timelike observers, the minimum will be positive. Cross-section is plotted for y = 0.

Fig. 20 — Eulerian stress-energy tensor components for the MTM. The T00 term is shown in Figure 18. Cross-section is plotted for y = 0.

Fig. 21 — Expansion and shear scalars for the MTM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

D. Lentz Metric (LM)

The Lentz metric (LM) is a radical departure from "top-hat" like spherically symmetric solutions proposed by Alcubierre and expanded on by Van Den Broeck, Bobrick, and Martire. The Lentz solution proposes to avoid WEC violations entirely, even in superluminal regimes. This is accomplished by recognizing that one issue with the prior warp metrics is the use of a single shift vector component, which provides a source of the negative energy. Lentz proposes that additional shift components can be used to construct a "potential" in the form of a linear wave equation: in Eq. (36), the sum of the second derivatives of the potential φ along x and y, minus the second derivative along z scaled by the inverse square of the characteristic velocity, equals the source density ρ,

which is related to the shift vector as

Nᵢ = ∂ᵢ φ (37)

and related to the Eulerian energy density by Eq. (38), in which the energy density E is one sixteenth of π times a combination of the second derivative of φ along z multiplied by the source density, plus the square of that derivative scaled by the inverse square of the characteristic velocity, minus four times the square of the mixed second derivative of φ along z and x.

A proposed solution to Equation 38 is a rhomboid-like shift vector along any two dimensions (in Cartesian coordinates). The APL team has built a version of this metric and evaluated it using the Warp Factory. Our form of the LM shift vector is shown in Figure 22. The APL evaluation of the energy density, shown in Figure 23, shows a similar result from Lentz's original work, with no major points of the energy density outside of the corners. This result differs from Lentz in that negative energy does exist for a few points, in addition to positive energy. However, this is likely due to the smoothing used in the shift vectors around the edges. (The smoothing applied by Lentz is unclear from his paper, but it is likely the case that the negative energy at the corners can be removed by more clever management of the gradients.) While the Eulerian energy density is mostly zero, the evaluations of the NEC and WEC show that violations occur throughout the structure, shown in Figure 24. The reason for this violation can be seen by taking a closer look at the stress-energy tensor, shown in Figure 25.

The culprit for WEC violation is the existence of Eulerian pressure terms at all shift vector boundaries with an absence of energy density. When transformed into other timelike observer frames, this appears as negative energy density to that observer. Finally, the scalars for the LM are shown in Figure 26.

Lentz's proposed soliton solution does provide a unique metric form that differs greatly from the "top-hat" form of Alcubierre and removes the major violations of Eulerian energy density. However, when running a full analysis across different timelike observers, this solution still results in WEC violations.

Fig. 22 — Metric functions for the LM. Only the shift vectors along X and Z have values different from Minkowski. These functions are shown in the comoving frame. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 23 — Eulerian energy density for the LM. Each point of energy density aligns with the corners of the rhomboid shape of the shift vector. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 24 — NEC and WEC for the LM. Direction of motion is along +Z. Plotted are the minimum values across all observers. When only positive energy density is seen across all null or timelike observers, the minimum positive value is shown. Cross-section is plotted for y = 0.

Fig. 25 — The non-zero stress-energy tensor components, in the Eulerian frame, for the LM. The T00 component is shown in Figure 23. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 26 — Expansion and shear scalars for the LM. Direction of motion is along +Z. Cross-section is plotted for y = 0.

E. Summary

The prior results of evaluating common and newly proposed metrics using Warp Factory are summarized in this section in terms of the statements of physicality (see Section III.C). The results for each metric are shown in Table 2. Using Warp Factory, elements of the stress-energy tensor can be fully explored to understand better how NEC and WEC violations occur.

Table 2 — A summary of the physicality conditions of the warp drive metrics discussed and evaluated in this paper. For the LM, the smoothing applied in this work does see regions of negative Eulerian energy, but this is due to not being able to replicate the smoothing applied in Lentz's paper. For the VDM, the reasonable mass condition could possibly be met, but requires very small regions with very large expansions which have not yet been analyzed by our code.

| Conditions | Alcubierre Metric | Van Den Broeck Metric | Bobrick-Martire Modified Time | Lentz Metric | |---|---|---|---|---| | Non-negative Eulerian energy | ✕ | ✕ | ✕ | ? | | Satisfies NEC everywhere | ✕ | ✕ | ✕ | ✕ | | Satisfies WEC everywhere | ✕ | ✕ | ✕ | ✕ | | Uses reasonable mass | ✕ | ? | ✕ | ✕ |

Across all of the metrics discussed in this paper, the common source of energy condition violation appears to stem from either the existence of negative energy density (as evaluated by a typical Eulerian observer) or the existence of pressure and momentum density terms exceeding the energy density at certain points. Once transformed into other timelike observers, these other terms become a negative energy seen by these other observers.

This leads to an important insight for warp analysis. To properly check for physicality, a full sampling of the null and timelike observers must be done to establish if any violation occurs. Determining only an Eulerian observer's energy density as non-negative is not a sufficient condition when discussing physicality.

VI. Optimization Results

The optimization method, discussed in Section IV.C, can take a starting metric and perturb it to generate an entirely new metric. This is a computationally intensive process but does offer some insight into what changes to a metric tend to improve physicality. To demonstrate the optimization capability a starting AM is used with a passenger volume radius of 20 m and a velocity in the z direction of 0.5c, with the bubble parameter σ set to 0.1. A genetic algorithm is run, with parameters described in Table 3. These optimizations are performed on an axisymmetric 2D coordinate system and translated following the optimization into 3D for visualization.

Table 3 — Optimization parameters run on the starting AM.

| Parameter | Value | Description | |---|---|---| | Passes | 10 | Number of passes that a given perturbation is allowed to add if an improvement is made. | | Generations | 540 | Number of optimization runs over a population. | | Population | 18 | Number of parallel optimizations run during a given generation. | | Perturbation Terms (6+4) | α, β2, β3, γ2, γ3, κ2 | Allowed metric functions to vary. |

The fitness metric used is the NEC, which takes all negative values of ν_NEC and integrates them across the entire volume. The starting metric has a fitness value of 7.2 × 10⁴⁴. (This is negative in nature, but the absolute value is used, with any reduction of this value being the improvement factor.) Following the optimization run, the fitness is improved by a factor of three to 2.1 × 10⁴⁴. Since the optimization can vary lapse, shift vectors, expansion, and shear terms, the final metric generated differs significantly from the starting AM. The resulting metric functions are shown in Figure 27. The Eulerian energy density for the starting and optimized metric is shown in Figure 28 and 29 respectively. The key difference is a contraction of the bubble energy density which is also accompanied by the addition of positive energy regions beyond it. A similar impact is seen in the energy condition violations as they are also pushed into the center region, shown in Figures 30 and 31.

Fig. 27 — Metric functions after optimization. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Looking at the full stress-energy components, shown in Figure 32, a much more intricate structure emerges with complex momentum flows and pressure regions. This complexity is partly due to the nature of the Gaussian perturbations themselves, which tend to cause "noisy" patches. A Gaussian function can approximate any kind of shape, but with a Monte Carlo approach, it can take a very long time to smooth out the wrinkles toward the final optimum solution, which itself is hard to determine in such a complex phase space.

Finally, comparing the scalars in Figures 33 and 34, the optimization resulted in the addition of concentric regions of the expansion and contraction. The inner region maintains the general expansion behind and contraction in front being slightly flattened along +X, but an outer region emerges which is opposite in sign. Again, it should be noted that the existence of expansion and contraction is not a requirement for a warp drive but can exist depending on the nature of the metric. Also differing from the AM is the shear, which again is everywhere positive but sees an addition of an inner and outer ring.

Fig. 28 — Eulerian energy density of the starting metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 29 — Eulerian energy density of the optimized metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 30 — NEC and WEC for the starting metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 31 — NEC and WEC for the optimized metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 32 — Eulerian stress-energy tensor components for the optimized metric. The T00 component is shown in Figure 29. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 33 — Expansion and shear scalars for the starting metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

Fig. 34 — Expansion and shear scalars for the optimized metric. Direction of motion is along +Z. Cross-section is plotted for y = 0.

The optimization performed here is only an example case to demonstrate its function. The method can be extended through the use of different fitness functions. For example, the NEC is useful in removing negative energy density, but often it results in the addition of positive energy to the world. Furthermore, the world and passenger volume boundaries restrict the possibility to remove global violations. If the defined world inherently has violations at the boundary (towards infinity), then the optimization is unable to change the root cause of the violation, often adding undesirable changes, such as additions of positive matter up to the world boundary. In general, the optimizer is best at removing local violations or optimizing the energy needed. Care must be taken to ensure that trivial solutions are not being generated. Certainly, a metric that does solve these issues at the onset can likely be optimized for its required energy, which will be an important step in addressing the reasonable mass requirement of a practical warp drive.

VII. Conclusion

Development of the Warp Factory codebase offers a powerful tool for the development of warp drives. The inherent challenges of solving the EFE present a hurdle to exploring novel warp metrics. The primary issue in developing warp metrics is overcoming the problems of its physicality, which requires a complicated analysis of timelike observers on the energy tensor components. To date, most analytical work considers a limited case of observers, which, for complicated metrics, might hide certain violations. The physicality of warp drives is not just a requirement of positive energy density, it is also the relation between momentum density, pressures, and shears within the stress-energy tensor. Evaluating this in a numerical framework allows users to establish and understand the dynamics of various warp metrics in a fully general manner. Warp Factory provides a tool for the visualization and analysis of arbitrary metrics with a focus on the evaluation of energy conditions, which is a unique capability when compared to most numerical GR tools.

In this paper, the functionality of Warp Factory has been demonstrated through the evaluation of classic and recent warp metrics, providing unique 3D analysis and visualizations of their stress-energy tensors, energy conditions, and scalars. One important insight is that the Lentz metric, which has been proposed as possibly solving the WEC violation, does still violate the WEC when analyzed across the space of timelike observers. Investigations into the impact of different passenger volumes on energy requirements and energy condition violations found that the overall shape has a marginal impact.

Warp Factory also implements an optimization module that provides a unique approach to warp research. While the method is computationally intensive and prone to limitations due to the large parameter space being explored, it offers insight into what elements of the warp metric tend to improve the energy conditions. It also allows changes to all metric functions as decided by the user, creating complicated structures not easily found or expressed analytically. The initial application of this to the AM sees some improvement in reducing violations through unique changes to the metric functions.

The APL team's next steps for Warp Factory are to expand the energy condition evaluation to include dominant and strong conditions while preparing the codebase for public use on Github. In addition, the team is developing functions for geodesic solvers to compute light rays and matter trajectories around a warp drive in a general 3D case. The team anticipates publishing follow-up papers on the results of these efforts.

Acknowledgments

Special thanks to the entire APL team and Sarah Dangelo for review and comments.

(The ten-item reference list is omitted for length; the complete text is at the source.)