A robophysical model of spacetime dynamics

Abstract

Systems consisting of spheres rolling on elastic membranes have been used to introduce a core conceptual idea of General Relativity: how curvature guides the movement of matter. However, such schemes cannot accurately represent relativistic dynamics in the laboratory because of the dominance of dissipation and external gravitational fields. Here we demonstrate that an “active” object (a wheeled robot), which moves in a straight line on level ground and can alter its speed depending on the curvature of the deformable terrain it moves on, can exactly capture dynamics in curved relativistic spacetimes. Via the systematic study of the robot’s dynamics in the radial and orbital directions, we develop a mapping of the emergent trajectories of a wheeled vehicle on a spandex membrane to the motion in a curved spacetime. Our mapping demonstrates how the driven robot’s dynamics mix space and time in a metric, and shows how active particles do not necessarily follow geodesics in the real space but instead follow geodesics in a fiducial spacetime. The mapping further reveals how parameters such as the membrane elasticity and instantaneous speed allow the programming of a desired spacetime, such as the Schwarzschild metric near a non-rotating blackhole. Our mapping and framework facilitate creation of a robophysical analog to a general relativistic system in the laboratory at low cost that can provide insights into active matter in deformable environments and robot exploration in complex landscapes.

Introduction

Systems consisting of spheres rolling on curved surfaces^1,2 are a well-known non-hydrodynamic analog to gravity. In such readily accessible systems, researchers have made intriguing connections to gravity such as Kepler-like laws, precession, and the stability of orbits. However, their studies have also found that these systems do not exactly mimic astrophysical gravity. For instance, the scaling between the period and radius is $T \propto r^{2/3}$³ instead of $T \propto r^{3/2}$ in Kepler’s third law. Additionally, the sphere on the elastic membrane is passive; as a result, not only do trajectories decay, but also the tunable parameters are limited to only the boundary conditions and the mass of the sphere.

Figure 1

A passive object (a marble) versus an active object (a robot) on a deformable membrane. While a passive marble is only subjected to friction and Earth’s gravity that leads to energy dissipation as $\vec {F}\cdot \vec {v}\propto \vec {a}\cdot \vec {v}<0$, an active object with an additional drive force can maintain steady-state motion with prescribed speed as $\vec {a}\cdot \vec {v}=0$.

We hypothesized that making the object “active” – an internally driven robot – would allow mechanical systems to better model GR in part because of the ability to study steady states. We further reasoned that the programmability and sensory capabilities of increasingly low-cost and powerful “robophysical” models^4,5 could allow the tuning of parameters that lead to inexact mimics of GR in passive systems. Indeed, our recent work⁶ built a framework to understand how the field-mediated dynamics of active agents on flexible membranes demonstrate in the words of Wheeler’s famous aphorism: “Matter tells spacetime how to curve and spacetime tells matter how to move”⁷. In particular, we showed that the spacetime followed by an active object (a wheeled circular robot) can be tuned by varying system parameters such as the membrane elasticity and the speed of the object.

Here we amplify on and extend the scheme introduced in⁶ and demonstrate how the activity can lead to an exact mapping to GR. We first show how an active object with a prescribed speed on an elastic membrane produces longer and more controllable trajectories compared with a passive marble. We then deduce the spacetime it follows, and subsequently program the spacetime with a Schwarzschild orbit as an example. We posit that a robot controlled in the way we describe could mechanically mimic blackhole dynamics in the laboratory at a low cost.

An active object with fixed speed on an elastic membrane

We first consider an active object prescribed with a constant speed on a circular elastic membrane. Later, we will discuss the general case of time-varying speed. To prevent the object from simply following a near-straight-line spatial geodesic with a spatial curvature

$$\begin{aligned} ds^2=\Psi ^2 dr^2+r^2 d\varphi ^2 \end{aligned}$$

(1)

where $\Psi ^2=1+z'^2$ and prime denotes the derivative with respect to r, the object must turn according to the instantaneous local slope, $-\nabla z$, the negative gradient of the terrain (membrane) height z. We enable this behavior by using a vehicle with a differential drive, a mechanism that permits a difference between the speeds of the vehicle’s two wheels to respond to the terrain slope while maintaining a center-of-mass speed prescribed by the motor. The vehicle drives straight on a leveled flat ground. When driving on a tilted flat surface with a constant gradient everywhere, the vehicle turns to align with this constant slope. In a more general case where the vehicle drives on a solid terrain with spatially varying gradient, it responds instantaneously by aligning with the local gradient. Finally, in our setup, the vehicle responds to the local gradient while the terrain (membrane) updates its shape with the evolving position of the vehicle, mimicking the interaction between matter and spacetime⁷. The membrane shape is affected by three factors. Earth’s gravity causes the membrane to sag with a parabolic profile due to the weight of the membrane itself. The tautness from boundary conditions such as the depth of the central depression of the membrane competes with the sagging. Finally, the local deformation from the vehicle creates an additional dimple in the membrane. We note that the sagging from the membrane’s weight and the deformation from the vehicle are analogous to the bending of background spacetime and the bending due to the moving object respectively in the context of GR.

Our terrain consists of a spandex membrane fixed to a circular frame with a diameter of 2.4 m. The central depth of the membrane is controlled by an actuator that adjusts the height of the membrane center. The cylindrical chassis of the vehicle is 3D printed and its instantaneous position is tracked by an overhead camera systems (Optitrack).

Figure 2

Trajectories of active and passive objects on an elastic membrane. (a) Sample perspective views of an active vehicle and a passive marble moving on a Spandex membrane. The time interval between two consecutive snapshots is 0.17 s. (b) The experimental trajectories, radius evolution, and speed evolution of the active (red) and the passive (blue) objects with the same mass (150 g) started from the same initial position and velocity on the same membrane. See the supplementary movie S1 for videos. (c) The simulation counterparts of (b).

We first compare the trajectories of the active vehicle with those of a passive marble having the same mass as the vehicle. We released the vehicle by placing it on the membrane after being turned on, and released the marble by placing it at the start of a guiding track. The velocity of the vehicle and marble were kept the same upon release by adjusting the voltage on the motor and the releasing height on the guiding track (Fig. 2a), respectively. The trajectories collected from experiments showed that the active vehicle produced much more persistent trajectories (Fig. 2c) than the passive marble, which barely completed one revolution (Fig. 2b).

To understand these orbits, we follow the models in^6,8. While a passive marble dissipates energy as $\vec {a}\cdot \vec {v}<0$ (Fig. 1), an active object can conserve its speed when the driving force dynamically balances with the friction and exactly makes $\vec {a}\cdot \vec {v}=0$ (Fig. 3a). Therefore, the acceleration for a constant-speed motion can be written as

$$\begin{aligned} \frac{a_{\varphi }}{r}= & {} \ddot{\varphi } + \frac{2\,\dot{r} \,{\dot{\varphi }}}{r} = \frac{a}{r}\cos {\theta } \end{aligned}$$

(2)

$$\begin{aligned} a_r= & {} \ddot{r} - \frac{r \,\dot{\varphi }^2}{\Psi ^2} + \frac{\Psi '}{\Psi }\dot{r}^2 = -\frac{a}{\Psi }\sin {\theta }, \end{aligned}$$

(3)

where $\theta$ is the heading angle between the radial direction and the velocity on an isotropic circular membrane.

Though the speed is constant, the change of the velocity (the scalar acceleration a) depends on the local slope $\gamma$ (Fig. 1). Since $\gamma$ varies with radius (position) r, a is also a function of r. Additionally, a should depend on velocity in general. However, given that the velocity has constant magnitude as the speed is constant, this dependence is reduced to one degree of freedom. For our convenience, we chose the direction of the velocity, $\theta$. If we consider an active object without chiral bias such that its trajectory has a mirror symmetry, the dependence of a (thus $a_{\varphi }$) on $\theta$ should be anti-symmetric about $\theta =0$, as otherwise the clockwise ($\theta (t=0)=\theta _0$) and counterclockwise ($\theta (t=0)=-\theta _0$) trajectories (Fig. 3b) will not be mirror reflections with each other. A first-order approximation with this symmetry could be $a\propto k(r)\sin {\theta }$ where the k(r) is the radial dependence due to the local slope $\gamma (r)$ that changes with radius. One could imagine k increases with the local slope $\gamma$. The detailed relation between k and $\gamma$ depends on the mechanical structure of the active object, but one can always Taylor expand this dependence. For preliminary study, here we assume linear dependence $k=C\gamma$.

Figure 3

Dynamics of the active vehicle. (a) The acceleration of an active vehicle is perpendicular to its velocity $\vec {v}$. (b) A non-chiral vehicle with a mirror-reflected initial velocity ${\mathcal {P}}\vec {v}$ will produce a mirror-reflected trajectory.

While an active object follows Eqs. (2, 3), a passive marble rolling on the membrane without slipping has a Lagrangian⁸

$$\begin{aligned} {\mathcal {L}}=\frac{7}{10}m\left( (1+z'(r)^2){\dot{r}}^2+r^2\dot{\varphi }^2\right) -mgz(r). \end{aligned}$$

(4)

We consider the Coulomb rolling friction $-\mu mg \hat{\varvec{v}}$ for a more realistic model. We neglect the higher-order effect ($O(\gamma ^2)$) from the slope on the normal force for simplicity. Plugging the corresponding dissipation function⁹ $D=\mu mg v$ where $v=\sqrt{{\dot{r}}^2+r^2\dot{\varphi }^2}$ into the Euler-Lagrange equation $\frac{d}{dt}\left( \frac{\partial {\mathcal {L}}}{\partial {\dot{q}}}\right) -\frac{\partial {\mathcal {L}}}{\partial q}=-\frac{\partial D}{\partial {\dot{q}}}$ where q is r or $\varphi$, we arrive at

$$\begin{aligned} \left( 1+z'^2\right) \ddot{r}+z'z''{\dot{r}}^2-r\dot{\varphi }^2+\frac{5}{7}gz'= & {} -\frac{5}{7}\mu g\frac{{\dot{r}}}{v} \end{aligned}$$

(5)

$$\begin{aligned} r^2\ddot{\varphi }+2r{\dot{r}}\dot{\varphi }= & {} -\frac{5}{7}\mu g\frac{r^2\dot{\varphi }}{v} \end{aligned}$$

(6)

The left-hand sides of the above equations are the same as the dynamical equations in¹ while the right hand sides correspond to the friction force.

Integration of the above models for the active vehicle and passive marble (Fig. 2c) shows qualitative agreement with the experiments. Figure 2c shows the integration of the active dynamics Eqs. (2, 3) and the passive dynamics Eqs. (5, 6) on the same simulated membrane with parameters measured from experiments. The active and passive objects started with the same position and velocity. The physical parameters are measured from experiments. The acceleration dependence on radius k for the active vehicle uses $k=C\gamma =C\partial _r z$ where z(r) is measured from the height of the static vehicle placed at different radii r. The proportionality C uses the ratio between acceleration and the gradient $\partial _r z$ at the radius close to the edge of the elastic membrane. We probe the friction coefficient for the passive marble by measuring the dissipation of mechanical energy in a designed experiment (see section “Probing the effective friction” of the Methods section).

Finding the spacetime corresponding to the orbits

A functional understanding of the orbital features (e.g., period, precession rate) would facilitate the creation of an intelligently programmed vehicle. One tool we can use to obtain such insights is the spacetime metric that describes these orbits. Any similarities between the experimentally inferred metric and known metrics could prove useful in understanding how the orbital features depend on system parameters.

The orbital dynamics we wish to map could be described by a diversity of metrics. We propose a simple but general metric that is analogous to GR in the weak-field limit and encodes the axi-symmetry of the system. Our metric has the form

$$\begin{aligned} ds^2 = -\alpha ^2 dt^2 +\Phi ^2\left( \Psi ^2 dr^2 + r^2 d\varphi ^2\right) \end{aligned}$$

(7)

with $\alpha = \alpha (r)$, $\Phi = \Phi (r), \Psi ^2=1+z'^2$. Here, the elements of the metric $g_{\alpha \beta }$ are zero except $g_{tt} = -\alpha ^2$, $g_{rr} = \Phi ^2\Psi ^2$ and $g_{\varphi \varphi } = \Phi ^2r^2$. Inserting $g_{\alpha \beta }$ into the Christoffel symbols $\Gamma ^{\mu }_{\alpha \beta }$ in the geodesic equations ${\mathop {x}\limits ^{\circ \circ }}^{\mu }+\Gamma ^{\mu }_{\alpha \beta }\mathring{x}_{\alpha }\mathring{x}_{\beta }=0$, we arrive at

$$\begin{aligned}{} & {} {\mathop {t}\limits ^{\circ \circ }} + \frac{\left( \alpha ^2\right) '}{\alpha ^2}\mathring{t}\mathring{r}= \frac{1}{\alpha ^2}\left( \alpha ^2\mathring{t}\right) ^\circ = 0 \end{aligned}$$

(8)

$$\begin{aligned}{} & {} {\mathop {\varphi }\limits ^{\circ \circ }} + \frac{\left( \Phi ^2\,r^2\right) '}{\Phi ^2\,r^2}\mathring{\varphi }\mathring{r}= \frac{1}{\Phi ^2\,r^2}\left( \Phi ^2\,r^2\mathring{\varphi }\right) ^\circ =0 \end{aligned}$$

(9)

$$\begin{aligned}{} & {} {\mathop {r}\limits ^{\circ \circ }} + \frac{(\alpha ^2)'}{2\Phi ^2\Psi ^2}\mathring{t}^2+ \frac{\left( \Phi ^2\Psi ^2\right) '}{2\Phi ^2\Psi ^2}\mathring{r}^2 -\frac{\left( \Phi ^2r^2\right) '}{2\Phi ^2\Psi ^2}\mathring{\varphi }^2 = 0 \end{aligned}$$

(10)

with $\lambda$ as an affine parameter and $\mathring{q}=dq/d\lambda ,{\mathop {q}\limits ^{\circ \circ }}=d^2q/d\lambda ^2$. From Eqs. (8, 9), we have that

$$\begin{aligned} \alpha ^2\mathring{t} = E = \textrm{constant}, \end{aligned}$$

(11)

$$\begin{aligned} \Phi ^2 r^2\mathring{\varphi } = L = \textrm{constant}. \end{aligned}$$

(12)

Both are consequences of conservation. Eq. (11) describes the conservation of energy, while Eq. (12) describes the conservation of angular momentum.

Using $\mathring{q}=(dq/dt)(dt/d\lambda )=\mathring{t}{\dot{q}}$ (see setion. “Converting derivatives” of the Methods section for details), the geodesic equations can be rewritten as

$$\begin{aligned} \ddot{\varphi } + \frac{2\dot{r} {\dot{\varphi }}}{r}= & {} \left[ \frac{(\alpha ^2)'}{\alpha ^2} - \frac{(\Phi ^2)'}{\Phi ^2}\right] \dot{r}\,{\dot{\varphi }} \end{aligned}$$

(13)

$${\ddot{r}} - \frac{{r\dot{\varphi }^{2} }}{{\Psi ^{2} }} + \frac{{\Psi ^{\prime } }}{\Psi }{\dot{r}^{2}} = \left[ {\frac{{\left( {\alpha ^{2} } \right)^{\prime } }}{{\alpha ^{2} }} - \frac{{\left( {\Phi ^{2} } \right)^{\prime } }}{{\Phi ^{2} }}} \right]\dot{r}^{2} + \frac{1}{{2\Phi ^{2} \Psi ^{2} }}\left[ {\left( {\Phi ^{2} } \right)^{\prime } v^{2} - \left( {\alpha ^{2} } \right)^{\prime } } \right]$$

(14)

where primes denote differentiation with respect to r.

Notice that the left-hand side of Eqs. (13, 14) are the components of the acceleration, $a_{\varphi }$ and $a_r$ respectively, from Eqs. (2, 3). When we substitute $\cos {\theta }={\dot{r}}/v,\sin {\theta }=r\dot{\varphi }/v$ and $a=k\sin {\theta }$ into Eqs. (2, 3), we find

$$\begin{aligned} \ddot{\varphi } + \frac{2\,\dot{r} \,{\dot{\varphi }}}{r}= & {} \frac{k}{v^2}{\dot{r}}\dot{\varphi } \end{aligned}$$

(15)

$$\ddot{r} - \frac{r \,\dot{\varphi }^2}{\Psi ^2} + \frac{\Psi^{\prime}}{\Psi }\dot{r}^2 = -\frac{k}{\Psi }\frac{r^2\dot{\varphi }^2}{v^2}.$$

(16)

Thus, comparing the right-hand sides of Eqs. (13, 14) and Eqs. (2, 3) and noticing that ${\dot{r}}^2+r^2\dot{\varphi }^2=v^2$ in Eq. (16), we obtain the following relationships between the metric functions $\alpha$ and $\Phi$ in terms of the speed of the vehicle and k.

$$\begin{aligned} \frac{(\alpha ^2)'}{\alpha ^2}= & {} \frac{k\Psi }{v^2}\left[ \frac{\Phi ^2v^2}{\alpha ^2-\Phi ^2v^2}\right] \end{aligned}$$

(17)

$$\begin{aligned} \frac{(\Phi ^2)'}{\Phi ^2}= & {} \frac{k\Psi }{v^2}\left[ \frac{2\Phi ^2v^2-\alpha ^2}{\alpha ^2-\Phi ^2v^2}\right] . \end{aligned}$$

(18)

Integration of the above equations yields

$$\begin{aligned} \alpha ^2= & {} -\frac{1}{C_1 v^2} + C_2 \cdot e^{-K/v^2} \end{aligned}$$

(19)

$$\begin{aligned} \Phi ^2= & {} \frac{\alpha ^2}{v^2}+C_1 (\alpha ^2)^2 \end{aligned}$$

(20)

where $K=K(r) \equiv \int _0^r k(s) \Psi (s) ds$. To determine the constants, we make use of the normalization condition and the fact that the metric should be flat at $k\rightarrow 0$.

The metric (Eq. 7) gives us the normalization condition $-1 = -\alpha ^2 \mathring{t}^2 + \Phi ^2(\Psi ^2\mathring{r}^2 + r^2 \mathring{\varphi }^2)$. To exploit this condition, we must eliminate the $d/d\lambda$ in $\mathring{r}$ like Eqs. (11, 12). Using $\mathring{q}=(dq/dt)(dt/d\lambda )=\mathring{t}{\dot{q}}$, Eqs. (11, 12), and the fact that $v^2=r^2\dot{\varphi }^2+{\dot{r}}^2$, we have

$$\begin{aligned} \mathring{r}^2= & {} \left( \frac{E}{\alpha ^2}{\dot{r}}\right) ^2\nonumber \\= & {} \frac{E^2}{(\alpha ^2)^2}\frac{1}{\Psi ^2}(v^2-r^2\dot{\varphi }^2)\nonumber \\= & {} \frac{E^2}{(\alpha ^2)^2}\frac{1}{\Psi ^2}\left[ v^2-r^2\left( \frac{\alpha ^2}{E}\mathring{\varphi }\right) ^2\right] \nonumber \\= & {} \frac{E^2}{(\alpha ^2)^2}\frac{1}{\Psi ^2}\left[ v^2-r^2\left( \frac{\alpha ^2}{E}\frac{L}{\Phi ^2 r^2}\right) ^2\right] . \end{aligned}$$

(21)

Plug the $\mathring{t},\mathring{r},\mathring{\varphi }$ derived above into the normalization condition, we now have

$$\begin{aligned} -1= & {} -\frac{E^2}{\alpha ^2} + \frac{\Phi ^2 E^2 v^2}{(\alpha ^2)^2}. \end{aligned}$$

(22)

Plugging in the $\alpha ^2$ and $\Phi ^2$ derived earlier (Eqs. 19, 20), we have $-\frac{1}{E^2}=C_1 v^2$, and therefore

$$\begin{aligned} C_1=-\frac{1}{v^2 E^2} \end{aligned}$$

(23)

Now, as promised earlier, we further determine the constants by making the metric flat when $k=0$. In fact, $k(r)=0$ indicates $K(r)=\int _{s=0}^r k(s)\Psi (s) = 0$. We set the lower limit of the integral to zero, without loss of generality, since otherwise the arbitrary constant will be absorbed by $C_2$. This limit reduces the metric to $\alpha _0^2 =-\frac{1}{C_1 v^2} + C_2$ and therefore $\Phi _0^2= \frac{\alpha _0^2}{v^2}+C_1 (\alpha _0^2)^2$ where $\alpha _0\equiv \lim _{k\rightarrow 0}\alpha$ and $\Phi _0\equiv \lim _{k\rightarrow 0} \Phi$. To satisfy the flatness that $\alpha _0=\Phi _0$, we need $\alpha _0^2=\frac{\alpha _0^2}{v^2}+C_1 (\alpha _0^2)^2$. By using $\alpha _0^2=1/C_1 v^2 +C_2$ again, we have

$$\begin{aligned} C_1C_2= & {} 1. \end{aligned}$$

(24)

The conditions in Eqs. (23, 24) settle the previously undetermined coefficients in the metric (Eqs. 19, 20). We finally arrive at

$$\begin{aligned} \alpha ^2= & {} E^2\left( 1-v^2 e^{-K/v^2}\right) \end{aligned}$$

(25)

$$\begin{aligned} \Phi ^2= & {} E^2 e^{-K/v^2} \left( 1-v^2 e^{-K/v^2}\right) . \end{aligned}$$

(26)

The constants E and L required to fix the metric have actual physical meanings. E is the constant energy associated with the fact that the metric is time-independent. L is the constant angular momentum associated with the metric’s $\varphi$-symmetry.

Thus our formulation indeed reveals that the vehicle does not simply follow spatial geodesics of the membrane but instead follows geodesics in an emergent spacetime (Eqs. 25, 26) generated by the global curvature, the local curvature, the active dynamics, and the differential mechanism. The resultant dynamics can now be understood as those of a test particle in a new spacetime where the active feature of the real particle, such as a persistently controlled speed, generates a non-splittable effective spacetime for the test particle (i.e. $g_{tt}$ is not constant). In the language of the work by Price¹⁰, the effects of curvature are now not restricted to space¹¹. That is, in general, the metric function $g_{tt}$ could depend on both the coordinate time (t) as well as the spatial coordinates. For a static metric (i.e., the metric functions are independent of time), the spacetime becomes splittable when $g_{tt}$ does not depend on the spatial coordinates. This leads to only spatial curvature. It was argued in¹⁰ that the spatial curvature is different from the spacetime curvature as it is devoid of gravity, i.e., a free particle initially at rest will remain at rest.

The essential contribution from the active drive is the persistent response to the local curvature, here enabled by the controlled constant speed unseen in passive systems. In fact, when the response of the turning to the local slope vanishes at the limit $v\rightarrow \infty$ such that $\alpha ^2 = \Phi ^2 = E^2(1-v^2)$, the metric Eq. (7) with components Eqs. (25, 26) reduces to a splittable (and flat) spacetime Eq. (1). On the other hand, when v is finite and controllable, the active locomotion provides more flexibility and programmability in constructing the desired spacetime in GR. For instance, programming an active agent with acceleration magnitude k decreasing with radius r makes an orbit precess in the same direction as the orbit, while a k increasing with r makes an orbit precess in the opposite direction as the orbit⁶. Such flexibility and programmability are more challenging to implement in passive and dissipative agents studied in the previous works^1,8.

Programming an arbitrary spacetime with a speed-varying robot

The metric Eqs. (25, 26) has shown us how the parameters of the system change the spacetime and thus the orbit. Now we see how we can solve the inverse problem of programming the desired spacetime using the system parameters (e.g., k(r) and v(r)).

In metric Eqs. (25, 26), we can tune the speed and membrane elasticity to change the spacetime of the orbits. However, here the spatial and radial metric are not yet completely disentangled yet. To have two degrees of freedom such that we can indeed program the spacetime arbitrarily, one could introduce another degree of freedom. For instance, if we allow the speed v to vary with the radius r (physical instantiation could be achieved by inferring the radius from the instantaneous tilting angle $\gamma$), Eqs. (17, 18) with $\Psi ^2\approx 1$ give the requirement for the mapping as

$$\begin{aligned} \frac{(\alpha ^2)'}{\alpha ^2}-\frac{(\Phi ^2)'}{\Phi ^2}= & {} \frac{v'}{v}+\frac{k}{v^2} \end{aligned}$$

(27)

$$\begin{aligned} \frac{(\Phi ^2)'v^2-(\alpha ^2)'}{2\Phi ^2}= & {} -k. \end{aligned}$$

(28)

These two equations above give us the recipe to create a desired spacetime by changing the speed of the vehicle with radius. For a desired metric with spatial curvature $\Phi ^2(r)$ and temporal curvature $\alpha ^2(r)$, we can solve for the required membrane elasticity and object speed by plugging in the curvatures into these two equations. The solution (see section “Programming the metric” of the Methods section for details) is

$$\begin{aligned} v(r)^2=\left( \int _{r_1}^r f(r')\cdot \frac{(\alpha ^2)'(r')}{\Phi ^2(r')}dr'\right) /f(r) \end{aligned}$$

(29)

where

$$\begin{aligned} f(r)=-e^{\int _{r_1}^r -2\frac{(\alpha ^2)'(r')}{\alpha ^2(r')}+\frac{(\Phi ^2)'(r')}{\Phi ^2(r')}dr'}. \end{aligned}$$

(30)

For instance, if we plug in the Schwarzschild metric in isotropic coordinates where $\alpha ^2(r)=1-r_s/r,\Phi ^2(r)=(1-r_s/r)^{-1}$, we arrive at the required membrane elasticity k(r) and active object speed v(r) as shown in Fig. 4a. Analytically,

$$\begin{aligned} v(r)^2= & {} r_s\frac{(r-r_s)^2}{r^3}+C\left( \frac{r-r_s}{r}\right) ^3 \end{aligned}$$

(31)

$$\begin{aligned} k(r)= & {} \frac{r_s(r-r_s)(r+Cr+r_s-Cr_s)}{2r^4} \end{aligned}$$

(32)

where

$$\begin{aligned} C=\frac{v_0^2r_0^3}{(r_0-r_s)^3}-\frac{r_s}{r_0-r_s}. \end{aligned}$$

(33)

Here, $v_0$ is the vehicle speed at $r_0$ as the boundary condition. For instance, one can use the inner radius as $r_0$.

Figure 4

Creating orbits in Schwarzschild spacetime with a speed varying particle. (a) The speed and membrane elasticity’s dependence on radius to create a Schwarzschild blackhole with $r_s=3.1$ mm. The inset shows a precessing orbit with $A=0.3$ m using this prescription. (b) Precession angle $|\Delta \varphi _{\textrm{prec}}|$ as a function of inverse latus rectum. (c) The relation between the orbital period, T, and the semi-major-axis, A, follows Kepler’s third law as $T\propto A^{3/2}$. Insets in (b) and (c) show the trajectories around the data points. See the supplementary movieS1 for continuous evolution of the orbits.

Simulations using this prescription show features of Schwarzschild orbits such as the linear dependence of the precession angle in terms of the inverse latus rectum. For Schwarzschild orbits with small precession, the precession angle increases with the inverse latus rectum as $\Delta \varphi _{\textrm{prec}} = 6\pi G^2M/(c^2l)=3\pi r_s \ell ^{-1}$ where G is the gravitational constant, M is the mass of the star, c is the speed of light, and $\ell \equiv A~(1-e^2)$ is the latus rectum. We evaluate the semi-major-axis A and the eccentricity e using the minimum and maximum radii: $A=(r_{\textrm{max}}+r_{\textrm{min}})/2,~e=(r_{\textrm{max}}-r_{\textrm{min}})/(r_{\textrm{max}}+r_{\textrm{min}})$. Fig. 4b shows the precession angle $\Delta \varphi _{\textrm{prec}}$ as a function of the inverse of the latus rectum, $\ell ^{-1}$, from simulations given $v_0=v(r_0=0.05\textrm{m})=0.225~\mathrm {m/s}$ and $r_s=0.0031$ m. The curve qualitatively follows the linear relationship, with small deviation from the theory due to the large precession angle. By changing $(r_0,v_0)$, we can obtain larger angular momenta and thus larger orbits around the same blackhole. These orbits show a relation between period T and semi-major-axis A that follows Kepler’s third law (Fig. 4c).

To achieve this relation in experiments, a vehicle must actively vary its speed with radius and a membrane must have a radially-dependent elastic modulus. We can attach a tilt sensor to infer the radius and change the speed accordingly. To program the membrane with radially varying profile $k(r)=Cg|\partial _r z|$, here we consider a membrane with linear elasticity following the Poisson Equation $\nabla \cdot E \nabla z = P$ where P is the unit load from the membrane gravity. One possible way to obtain the desired k(r) is to create an elastic material with a radially varying thickness $P=P(r)$. Another option is to fabricate a membrane with a radially varying modulus $E=E(r)$.

Programming a spacetime with a non-diagonal metric

In principle, one can even create orbits from a non-diagonal metric by breaking symmetries. Here we show that breaking the axis symmetry can create a spacetime with nonzero $g_{t\phi }$, which is essential in perhaps the most well-known non-diagonal metric, the Kerr metric for a rotating blackhole. Experimentally, this could be done by adding tangential perturbation to the substrate. See Fig. 5 for an illustration.

Let us consider a $(2+1)$D metric $g_{\mu \nu }$ where the only off-diagonal term is the $g_{t\phi }$. Following the same methodology we used for the Schwarzschild metric, the geodesic equation is now

$$\begin{aligned} {\mathop {t}\limits ^{\circ \circ }}= & {} -\Gamma ^t_{tr}\mathring{t}\mathring{r}{{-\Gamma ^t_{r\phi }\mathring{r}\mathring{\phi }}}\end{aligned}$$

(34)

$$\begin{aligned} {\mathop {r}\limits ^{\circ \circ }}= & {} -\Gamma ^r_{tt}\mathring{t}^2{{-\Gamma ^r_{t\phi }\mathring{t}\mathring{\phi }}}-\Gamma ^r_{rr}\mathring{r}^2-\Gamma ^r_{\phi \phi }\mathring{\phi }^2 \end{aligned}$$

(35)

$$\begin{aligned} {\mathop {\phi }\limits ^{\circ \circ }}= & {{-\Gamma ^{\phi }_{tr}\mathring{t}\mathring{r}}}-\Gamma ^{\phi }_{r\phi }\mathring{r}\mathring{\phi } \end{aligned}$$

(36)

where $\Gamma_{r\phi}^t \mathring{r}\mathring{\phi},\Gamma_{t\phi}^r \mathring{t}\mathring{\phi},\Gamma_{tr}^{\phi}\mathring{t}\mathring{r}$ are from the off-diagonal components in addition to the diagonal metric (for instance, Schwarzschild) we worked on in the previous section. See section “Programming metric with off‑diagonal” of the Methods section for the technical details of this section.

The conserved quantities are now generalized to be $Q_t=P_{tt}(r)\mathring{t}+P_{t\phi }\mathring{\phi }$ and $Q_{\phi }=P_{\phi t}(r)\mathring{t} +P_{\phi \phi }\mathring{\phi }$ where

$$\begin{aligned} P^{-1}P'= \begin{pmatrix} \Gamma _{tr}^t&{\Gamma _{r\phi }^{t}}\\ {\Gamma _{{tr}}^{\phi }}&{}\Gamma _{r\phi }^{\phi } \end{pmatrix} \text {with } P\equiv \begin{pmatrix} P_{tt} &{} P_{t\phi }\\ P_{\phi t} &{} P_{\phi \phi } \end{pmatrix}. \end{aligned}$$

(37)

While $\Gamma _{r\phi }^t$ and $\Gamma _{tr}^{\phi }$ are now nonzero such that P is not simply $\text {diag}\{\alpha ^2(r),\Phi ^2(r)r^2\}$ as we have in the diagonal case, the conserved quantities $(Q_t,Q_{\phi })^T=P(\mathring{t},\mathring{\phi })^T$ still give us $(\mathring{t},\mathring{\phi })^T=P^{-1}(Q_t,Q_{\phi })^T\equiv (f(r),g(r))^T$ as functions of r. From this, we obtain the equations of motion with respect to time t instead of affine parameter $\lambda$ as

$$\begin{aligned} \ddot{\phi }+\frac{2{\dot{r}}\dot{\phi }}{r}= & {} c_{r\phi }{\dot{r}}\dot{\phi }-\Gamma _{tr}^{\phi }{\dot{r}} \end{aligned}$$

(38)

$$\begin{aligned} \ddot{r}-r\dot{\phi }^2= & {} c_{\phi \phi }\dot{\phi }^2-\Gamma _{t\phi }^r\dot{\phi }+c_0 \end{aligned}$$

(39)

where the c’s are all functions of r such that $c_{\phi \phi }=(f'/f+\Gamma _{rr}^r)r^2-\Gamma _{\phi \phi }^r-r$, $c_0=-\Gamma _{tt}^r-(f'/f+\Gamma _{rr}^r)v^2$, $c_{r\phi }=-\Gamma _{r\phi }^{\phi }-f'/f+2/r$.

To accommodate the new terms with $\dot{\phi }$ and ${\dot{r}}$, one could use a rotating object¹² (for instance a tilted slab shown in Fig. 5) to locally add an azimuthal perturbation $\delta$. This perturbation would break the symmetry such that the magnitude of acceleration $a=k (\sin {\theta } + \delta )$ has a bias in one chirality over the other since $a(\theta ')\ne a(-\theta ')$. Noting $\sin {\theta }=r\dot{\phi }/v$, the equation of motion of the vehicle on a substrate with broken axial symmetry is

$$\begin{aligned} \ddot{\phi }+\frac{2{\dot{r}}\dot{\phi }}{r}= & {} \frac{1}{r}\frac{{\dot{r}}}{v}k\left( \frac{r\dot{\phi }}{v}+\delta \right) \end{aligned}$$

(40)

$$\begin{aligned} \ddot{r}-r\dot{\phi }^2= & {} -\frac{r\dot{\phi }}{v}k\left( \frac{r\dot{\phi }}{v}+\delta \right) \end{aligned}$$

(41)

$\delta$ is a function of r, which can be determined by matching the r-dependent functions in Eqs. (38–39) and Eqs. (40–41) either analytically or numerically.

Figure 5

Proposed scheme to create a spacetime with a non-diagonal metric. A controlled deformation (e.g., a tilted slab) rotating about the central axis shown in (A) could plausibly induce an asymmetric dependence on heading $\theta$ for the acceleration magnitude in blue instead of the symmetric counterpart in red. The spacetime metric resulting from this proposed setup would have a non-diagonal component $g_{t\phi }$.

Discussion

In this work, we demonstrated how the use of an active particle – a wheeled robot – moving on an elastic membrane can generate a system which can mimic the dynamics of bodies in arbitrary spacetime. Given the flexibility in constructing and programming such robophysical devices, our system makes for an attractive target to push toward a mechanical analog GR system. While superficially our system resembles the educational tool used to motivate Einstein’s view of spacetime curvature influencing matter trajectories^1,3,8, unlike such systems which are not good analogs of GR, the activity allows the dynamics of the vehicle to be dictated by the curvature of “spacetime”, not just the curvature of space as in splittable spacetimes (where $g_{tt}$ is constant)¹⁰. As such, our system can be used as an experimental example of GR in upper-division physics courses to enhance students’ hands-on understanding^13,14 of orbits and curved spacetime. Thus we posit that mechanical analog “robophysical”^4,15 systems can complement existing fluid^16,17, condensed matter¹⁸, atomic, and optical^19,20,21 analog gravity systems²² given the ability to create infinite types of spacetimes. We might even generate analogies to wave-like systems^23,24,25; for example, one could increase the speed of the vehicle to be comparable to disturbance propagation (such that the membrane would follow the wave equation). There has been much theoretical progress on the analogy between spacetime and solid mechanics²⁶. We posit mechanical devices, such as the one we introduce in this work, could help explore the theoretical proposal with the reverse methodology. For instance, we can design systems to monitor the elasticity of known spacetimes and probe interesting problems such as the occurrence of topological defects and plastic deformations in spacetime.

We advocate for the potential research applications of our robotic vehicle as an analog system. In astrophysics, we typically are limited to an observational understanding of phenomena borne out of an inability to perform physical experiments with customized parameters. Though simulation provides an alternative, limitations such as sub-grid physics and runtime practicality for high-resolution or multiscale dynamics²⁷, path dependency and contingency of composite models²⁸, and the nature of modeling approximation when carrying out explorative research²⁹ are nontrivial obstacles that may hinder the validity and applicability of these simulations. Physical analog systems provide simple parameter variation, intellectual accessibility, and potentially offer ease of solution due to their classical nature²². With the flexibility of parameters introduced in this paper, the length and time scales not accessible in celestial systems are possible to access in their laboratory analogs. For instance, it is possible to probe physics around or even within the event horizon of a blackhole to probe theories³⁰ that are challenging to test in the original system. See section “Orbits around the horizon” of the Methods section for one such example of repulsive orbits inside the event horizon.

Beyond its role as a mechanical analog for GR, this framework could also provide a new perspective to understanding active matter undergoing field-mediated interactions^6,31. For instance, the spacetime metric of the agents’ motion can both guide our choice of parameter values to alter orbital features, like the precession direction, and influence our design of control schemes that accomplish tasks like helping multiple agents avoid mergers on the membrane⁶.

Methods

Probing the effective friction

We simplify the complicated rolling friction and membrane dissipation by employing an effective friction constant that absorbs all dissipative forces. We probe its magnitude by doing the following experiment. We release the marble at the rim of the circular membrane with zero speed and thus zero kinetic energy (Fig. 6). The marble then rolls radially towards the center, passes through the center, and stops before it reaches the other end of the diameter due to the effective rolling friction. Absorbing the loss of mechanical energy into the dissipation from the effective rolling friction $f_{\textrm{roll}}$ for a distance of $\ell$, we arrive at

$$\begin{aligned} f_{\textrm{roll}}\ell =mg\Delta h \end{aligned}$$

(42)

The measurements from experiment that $\ell = 1.5$ m, $\Delta h = 0.1$ m give the effective friction coefficient $\mu = f_{\textrm{roll}}/mg=\Delta h/\ell = 0.07 \sim 0.1$.

Figure 6

An experiment to probe the effective friction.

Converting derivatives

With the help of $\mathring{q}\equiv \frac{dq}{d\lambda } = \frac{dt}{d\lambda }\frac{dq}{dt}$ and $\alpha ^2\mathring{t} = E, \Phi ^2 r^2\mathring{\varphi } = L$ in Eqs. (11, 12), we have

$$\begin{aligned} \mathring{t}= & {} \frac{E}{\alpha ^2} \end{aligned}$$

(43)

$$\begin{aligned} {\mathop {t}\limits ^{\circ \circ }}= & {} \frac{dt}{d\lambda }\frac{d\mathring{t}}{dt}\nonumber \\= & {} \frac{E}{\alpha ^2} \frac{d}{dt}\left( \frac{E}{\alpha ^2}\right) \nonumber \\= & {} -\frac{E^2 (\alpha ^2)'}{(\alpha ^2)^3} {\dot{r}} \end{aligned}$$

(44)

$$\begin{aligned} \mathring{r}= & {} \frac{dt}{d\lambda }\frac{dr}{dt}=\frac{E}{\alpha ^2}{\dot{r}} \end{aligned}$$

(45)

$$\begin{aligned} {\mathop {r}\limits ^{\circ \circ }}= & {} \frac{dt}{d\lambda }\frac{d\mathring{r}}{dt}\nonumber \\= & {} \frac{E}{\alpha ^2}\cdot \frac{d}{dt}\left( \frac{E}{\alpha ^2}{\dot{r}}\right) \nonumber \\= & {} \frac{E^2}{(\alpha ^2)^2}\cdot \left( -\frac{(\alpha ^2)'}{\alpha ^2}{\dot{r}}^2+\ddot{r}\right) \end{aligned}$$

(46)

$$\begin{aligned} \mathring{\varphi }= & {} \frac{dt}{d\lambda }\frac{d\varphi }{dt}=\frac{E}{\alpha ^2}\dot{\varphi } \end{aligned}$$

(47)

$$\begin{aligned} {\mathop {\varphi }\limits ^{\circ \circ }}= & {} \frac{dt}{d\lambda }\frac{d\mathring{\varphi }}{dt}\nonumber \\= & {} \frac{E}{\alpha ^2 } \frac{d}{dt}\left( \frac{E}{\alpha ^2}\dot{\varphi }\right) \nonumber \\= & {} \frac{E^2}{(\alpha ^2)^2}\cdot \left( -\frac{(\alpha ^2)'}{\alpha ^2}{\dot{r}}\dot{\varphi }+\ddot{\varphi }\right) . \end{aligned}$$

(48)

Programming the metric

By eliminating the k in Eqs. (27, 28), we get

$$\begin{aligned} M V-V'=\frac{(\alpha ^2)'}{\Phi ^2} \end{aligned}$$

(49)

where $M(r)=2(\alpha ^2)'(r)/\alpha ^2(r)-(\Phi ^2)'(r)/\Phi ^2(r)$ and $V(r)=v^2(r)$.

We can multiply a function f(r) to both sides of Eq. (49) to make the left-hand side exact. Noting $(fV)'=f'V+fV'$, we need $f'/f=-M$. Therefore,

$$\begin{aligned} f(r)=-e^{\int _{r_1}^r -M(r')dr'}. \end{aligned}$$

(50)

With this f, we now have $(fV)'=f~(\alpha ^2)'/\Phi ^2$. So,

$$\begin{aligned} V(r)=\left( \int _{r_1}^r f(r')\cdot \frac{(\alpha ^2)'(r')}{\Phi ^2(r')}dr'\right) \cdot \frac{1}{f(r)}. \end{aligned}$$

(51)

By plugging in the Schwarzschild metric in isotropic coordinates $\alpha ^2(r)=1-r_s/r,\Phi ^2(r)=(1-r_s/r)^{-1}$, we have

$$\begin{aligned} f(r)= & {} -e^{\int _{r_1}^r -\frac{3r_s}{r'(r'-r_s)}dr'}\nonumber \\= & {} -e^{(C_1+3\log {(r/(r-r_s))})}\nonumber \\= & {} -C_2\cdot \left( \frac{r}{r-r_s}\right) ^3. \end{aligned}$$

(52)

Therefore,

$$\begin{aligned} V(r)= & {} \left( \int _{r_1}^r C_2\cdot \left( \frac{r'}{r'-r_s}\right) ^3\cdot \frac{(r'-r_s)r_s}{r'^3}dr'\right) /f(r)\nonumber \\= & {} \left( \int _{r_1}^r C_2\cdot \frac{r_s}{(r'-r_s)^2}dr'\right) /f(r)\nonumber \\= & {} \left( -C_2\frac{r_s}{r-r_s}+C_3\right) /\left( -C_2\cdot \left( \frac{r}{r-r_s}\right) ^3\right) \nonumber \\= & {} r_s\frac{(r-r_s)^2}{r^3}+C\left( \frac{r-r_s}{r}\right) ^3 \end{aligned}$$

(53)

$$\begin{aligned} k(r)= & {} -\frac{(\Phi ^2)'V-(\alpha ^2)'}{2\Phi ^2}\nonumber \\= & {} \frac{(r-r_s)r_s(r+Cr+r_s-Cr_s)}{2r^4}. \end{aligned}$$

(54)

To program the active object physically, we want to prescribe the speed $v_0$ at a certain radius (say the inner radius $r_0$) so that $V(r_0)=v_0^2$, we need

$$\begin{aligned} C=\frac{v_0^2r_0^3}{(r_0-r_s)^3}-\frac{r_s}{r_0-r_s}. \end{aligned}$$

(55)

Further, a reasonable speed $v_c$ at a characteristic orbit size (say the circular orbit $r_c$) will limit the size of the Schwarzschild radius $r_s$ (the size of the blackhole) with $\frac{V(r_c;r_s)}{r_c}=k(r_c)$.

Programming metric with off-diagonal

For metric with nonzero $g_{t\phi }$, the nonzero Christoffel symbols are

$$\begin{aligned} {\Gamma ^t_{tr}}= & {} {\frac{1}{2}(g^{tt}g_{tt,r}+g^{t\phi }g_{\phi t,r})} \end{aligned}$$

(56)

$$\begin{aligned} {\Gamma ^t_{r\phi }}= & {} {\frac{1}{2}(g^{tt}g_{t\phi ,r}+g^{t\phi }g_{\phi \phi ,r})} \end{aligned}$$

(57)

$$\begin{aligned} {\Gamma ^r_{tt}}= & {} {-\frac{1}{2}g^{rr}g_{tt,r}} \end{aligned}$$

(58)

$$\begin{aligned} {\Gamma ^r_{t\phi }}= & {} {-\frac{1}{2}g^{rr}g_{t\phi ,r}} \end{aligned}$$

(59)

$$\begin{aligned} {\Gamma ^r_{rr}}= & {} {\frac{1}{2}g^{rr}g_{rr,r}} \end{aligned}$$

(60)

$$\begin{aligned} {\Gamma ^r_{\phi \phi }}= & {} {-\frac{1}{2}g^{rr}g_{\phi \phi ,r}} \end{aligned}$$

(61)

$$\begin{aligned} {\Gamma ^{\phi }_{tr}}= & {} {\frac{1}{2}(g^{\phi \phi }g_{\phi t,r}+g^{\phi t}g_{tt,r})} \end{aligned}$$

(62)

$$\begin{aligned} {\Gamma ^{\phi }_{r\phi }}= & {} {\frac{1}{2}(g^{\phi t}g_{t\phi ,r}+g^{\phi \phi }g_{\phi \phi ,r})} \end{aligned}$$

(63)

The geodesic equation is therefore

$$\begin{aligned} {\mathop {t}\limits ^{\circ \circ }}= & {} -\Gamma ^t_{tr}\mathring{t}\mathring{r}-\Gamma ^t_{r\phi }\mathring{r}\mathring{\phi } \end{aligned}$$

(64)

$$\begin{aligned} {\mathop {r}\limits ^{\circ \circ }}= & {} -\Gamma ^r_{tt}\mathring{t}^2-\Gamma ^r_{t\phi }\mathring{t}\mathring{\phi }-\Gamma ^r_{rr}\mathring{r}^2-\Gamma ^r_{\phi \phi }\mathring{\phi }^2 \end{aligned}$$

(65)

$$\begin{aligned} {\mathop {\phi }\limits ^{\circ \circ }}= & {} -\Gamma ^{\phi }_{tr}\mathring{t}\mathring{r}-\Gamma ^{\phi }_{r\phi }\mathring{r}\mathring{\phi } \end{aligned}$$

(66)

Assume now the conserved quantities are $Q_t=P_{tt}\mathring{t}+P_{t\phi }\mathring{\phi },Q_{\phi }=P_{\phi t}\mathring{t}+P_{\phi \phi }\mathring{\phi }$, then

$$\begin{aligned} \frac{dQ_t}{d\lambda }=P'_{tt}\mathring{r}\mathring{t}+P_{tt}{\mathop {t}\limits ^{\circ \circ }}+P'_{t\phi }\mathring{r}\mathring{\phi }+P_{t\phi }{\mathop {\phi }\limits ^{\circ \circ }}= & {} 0 \end{aligned}$$

(67)

$$\begin{aligned} \frac{dQ_{\phi }}{d\lambda }=P'_{\phi t}\mathring{r}\mathring{t}+P_{\phi t}{\mathop {t}\limits ^{\circ \circ }}+P'_{\phi \phi }\mathring{r}\mathring{\phi }+P_{\phi \phi }{\mathop {\phi }\limits ^{\circ \circ }}= & {} 0 \end{aligned},$$

(68)

which can be also written as

$$\begin{aligned} P\begin{pmatrix} {\mathop {t}\limits ^{\circ \circ }}\\ {\mathop {\phi }\limits ^{\circ \circ }} \end{pmatrix} = -P'\begin{pmatrix} \mathring{r}\mathring{t}\\ \mathring{r}\mathring{\phi } \end{pmatrix}, \text {where } P\equiv \begin{pmatrix} P_{tt} &{} P_{t\phi }\\ P_{\phi t} &{} P_{\phi \phi } \end{pmatrix} \end{aligned}.$$

(69)

When

$$\begin{aligned} P^{-1}P'= \begin{pmatrix} \Gamma _{tr}^t&{}\Gamma _{r\phi }^t\\ \Gamma _{tr}^{\phi }&{}\Gamma _{r\phi }^{\phi } \end{pmatrix} \end{aligned}$$

(70)

is satisfied, Eq. (69) matches the t and $\phi$ components of the geodesic equations. When the metric is diagonal, it will reduce to the result for diagonal metric $P'/P=\text {diag}\{\Gamma _{tr}^t,\Gamma _{r\phi }^{\phi }\}$ since both P and $\Gamma$ are diagonal. The conserved quantities $P(\mathring{t},\mathring{\phi })^T=(Q_t,Q_{\phi })^T$ give us $\mathring{t}$ and $\mathring{\phi }$ as functions of r: $(\mathring{t},\mathring{\phi })^T=P^{-1}(Q_t,Q_{\phi })^T\equiv (f(r),g(r))^T$. This allows us to convert the geodesic equation in terms of affine parameter $\lambda$ to time t by using method shown in section “Converting derivatives”. Plug in the conversions $\mathring{r}=f{\dot{r}},{\mathop {r}\limits ^{\circ \circ }}=f(f'\dot{(}r)^2+f\ddot{r}),\mathring{\phi }=f\dot{\phi },{\mathop {\phi }\limits ^{\circ \circ }}=f(f'{\dot{r}}\dot{\phi }+f\ddot{\phi })$ into Eq. (65, 66), we arrive at

$$\begin{aligned} \ddot{\phi }+\frac{2{\dot{r}}\dot{\phi }}{r}= & {} c_{r\phi }{\dot{r}}\dot{\phi }-\Gamma _{tr}^{\phi }{\dot{r}} \end{aligned}$$

(71)

$$\begin{aligned} \ddot{r}-r\dot{\phi }^2= & {} c_{\phi \phi }\dot{\phi }^2-\Gamma _{t\phi }^r\dot{\phi }+c_0 \end{aligned}$$

(72)

where c’s are all functions of r that $c_{\phi \phi }=(f'/f+\Gamma _{rr}^r)r^2-\Gamma _{\phi \phi }^r-r$, $c_0=-\Gamma _{tt}^r-(f'/f+\Gamma _{rr}^r)v^2$, $c_{r\phi }=-\Gamma _{r\phi }^{\phi }-f'/f+2/r$.

Orbits around the horizon

We use the same model in previous section for Schwarzschild spacetime to integrate vehicle trajectory with controlled speed v(r) and acceleration response to designed membrane elasticity k(r). The integration uses

$$\begin{aligned} r\ddot{\phi }+2{\dot{r}}\dot{\phi }= & {} \left( \frac{v'}{v}+\frac{k}{v^2}\right) r{\dot{r}}\dot{\phi } \end{aligned}$$

(73)

$$\begin{aligned} \ddot{r}-r\dot{\phi }^2= & {} \left( \frac{v'}{v}+\frac{k}{v^2}\right) {\dot{r}}^2-k \end{aligned},$$

(74)

which is derived from

$$\begin{aligned} a_r= & {} a_t\cos {\theta }-a_n\sin {\theta } \end{aligned}$$

(75)

$$\begin{aligned} a_{\phi }= & {} a_t\sin {\theta }+a_n\cos {\theta } \end{aligned}$$

(76)

where the normal and tangential components of the acceleration are $a_n=k(r)\sin {\theta },a_t=dv/dt=(\partial v/\partial r) {\dot{r}}\equiv v' {\dot{r}}$. The orbits around and inside the event horizon are shown in Fig. 7.

Figure 7

Orbits around the horizon ($r_s=3.1$ mm). As the initial radius gets closer to the event horizon shown in dashed line, the orbit approaches the blackhole and eventually is captured by the singularity at $r=r_s$. Using the method described earlier to integrate orbits outside the horizon, the orbit inside the horizon starts to show intriguing repulsion as predicted by³⁰.