the time of flight is spent at lower radius, and closer to the Sun. Our intent
is not to perform an extensive analysis on the number of levels of sub-optimal
solutions, but these results indicate that formulation of minimum-fuel optimal
control problems using Cartesian coordinates will result, most of the time, in
sub-optimal solutions.

The existence of these sub-optimal solutions (associated with $N_{\text{rev}}$
en route revolutions) is a direct consequence of large time of flight and
thrust acceleration magnitude and the additional freedom to make as many
revolutions as possible, within the given time. Figures 66-68 show
trajectories and thrust profiles of the four sub-optimal solutions.

Figure 67: Earth-to-Dionysus sub-optimal solution with $m_{f}=2465$ kg.

Figure 68: Earth-to-Dionysus sub-optimal solution with $m_{f}=2672$ kg.

###  7.2 GTO To A Halo Orbit Around L1

This example illustrates how to design an optimal low-thrust transfer from a
Geostationary Transfer Orbit (GTO) to an L1 halo orbit in the restricted
three-body problem associated with the Earth-Moon system. The equations of
motion are given in [99]. The departure of the spacecraft is from the perigee
of an elliptical orbit with perigee $\times$ apogee altitude ratio of
$400\times 35,864$ km (GTO). The spacecraft with an initial mass of
$m_{0}=1500$ kg uses an electric engine with a specific impulse of
$I_{sp}=3000$ seconds.

Figure 69: Halo orbit and point of entrance to the orbit.

The target orbit is a halo orbit with an out-of-plane amplitude of 8000 km
[martin2010optimal]. Figure 69 shows the halo orbit and the point on orbit at
which the spacecraft will enter the orbit. The parameters used in numerical
simulations are : $g_{0}=9.80665$ (m/s2), Length unit (LU) = 3.84405000
$\times 10^{5}$ (km), Time unit (TU) = $3.75676967\times 10^{5}$ (s), Velocity
unit (VU) = $1.02323281$ (km/s), $\mu=1.21506683\times 10^{-2}$ (LU3/TU2). LU
in actually the reference semi-major axis of the Moon’s orbit relative to the
Earth (the distance between Earth and Moon).

Figure 70: Switching surface for sub-optimal trajectories for GTO-to-Halo
problem with $N_{\text{rev}}=6$.

This problem is more difficult compared to the previous cases where the
dynamics were simpler. For the three-body dynamical models, convergence is
found to be difficult to achieve using a finite difference approach for
calculating sensitivities. Therefore, the STM method (with hyperbolic tangent
smoothing) is used for the sensitivities according to the procedure described
in [71]. When integrating the state and co-state equations, we make use of
hyperbolic tangent smoothing to enhance the convergence through a continuation
procedure.

Finding the minimum-thrust solution is the first step for constructing the
switching surface. However, the minimum-thrust solution is not unique, which
is an important factor. For instance, Figures 72 and 72 show the switching
function and thrust profile and trajectory for a minimum-thrust solution with
$T_{\text{min}}\approx 9.675$ N with $N_{\text{rev}}=6$. Clearly, the
switching function is positive along the whole trajectory. However, this is
not the fundamental minimum-thrust solution since the number of revolutions
are greater than the global solution. If this local extremal solution is
considered as the base solution and the thrust magnitude is swept, Figure 70
shows the resulting switching surface for $T\in[9.675,20]$ N with
$N_{\text{rev}}=6$. Note that when Cartesian elements are used for formulating
the TPBVPs, the number of revolutions does not appear explicitly. There
appears to be seven main thrust ridges that remain at high thrust levels.
There is a late-departure boundary, whereas the last thrust ridge is attached
to the terminal time, which indicates that there is no early-arrival boundary.

Figure 71: GTO-to-Halo switching function with $T\approx 9.675$ N and
$N_{\text{rev}}=6$.

Figure 72: GTO-to-Halo trajectory with $T\approx 9.675$ N and
$N_{\text{rev}}=6$.

However, the global minimum-thrust solution of this problem has been found to
correspond to $T_{\text{min}}\approx 8.141$ N with $N_{\text{rev}}=5$. Figures
74 and 74 show the switching function and thrust profile and trajectory for
the fundamental minimum-thrust solution with $T_{\text{min}}\approx 8.141$ N.
It is interesting to note that in Figure 72, the profile of the switching
function has a different behavior during the time interval $t\in[0.5,1]$
compared to the smooth one in Figure 74.

Figure 73: GTO-to-Halo switching function with
$T\approx 8.141$ N.

Figure 74: GTO-to-Halo trajectory with $T\approx 8.141$ N.

Figure 76 shows the resulting optimal switching surface when $T\in[8.141,15]$
N is considered. There appears to be only seven main thrust ridges that remain
at “high” thrust levels. The boundary conditions result in a switching surface
without early-arrival or later-departure boundaries. Note that the thrust is
not necessarily high enough as the width of the thrust ridges are still large.

Figure 75: Switching surface and thrust profile for optimal trajectories,
$N_{\text{rev}}=5$ for GTO-to-Halo problem with $T\in[8.141,15]$ N.

Figure 76: Switching surface for optimal trajectories, $N_{\text{rev}}=5$ for
GTO-to-Halo problem $T\in[8.141,40]$ N.

Figure 76 shows the optimal switching surface for $T\in[8.141,40]$ N when
$\rho_{\text{min}}=9.15\times 10^{-6}$. Note that in practice, a relatively
small value is chosen for the smoothing parameter in order to facilitate the
generation of the switching surface. For particular regions of interest, it is
possible to use a smaller value for $\rho_{\text{min}}$.

Figure 77: GTO-to-Halo switching function and thrust profile with $T\approx
24.65$ N.

Figure 78: GTO-to-Halo trajectory with $T\approx 24.65$ N.

There is an interesting bifurcation phenomenon at $T\approx 19$ N, see Figure
76. The nature of this bifurcation point is completely different from the
previously shown bifurcations in the two-body dynamics (where a corner point
was the beginning point of the bifurcation and the beginning of the thrust arc
was attached to a thrust ridge, see Fig. 10). In addition, in the two-body
dynamics $\ddot{S}>0$ in the vicinity of the bifurcation point, whereas here
the sign of $\ddot{S}$ changes with $S=\dot{S}=\ddot{S}=0$ at the bifurcation;
hence, this bifurcation appears in the middle of a coast canyon. Clearly, the
creation of this thrust ridge has also visibly correlated to sharp features on
the two closest thrust ridges to its right and left where two corner points
are noticeable and their time duration (width) shrinking rate with increasing
thrust is increased.

Figures 78 and 78 show the switching function and the trajectory for this
bifurcation point. In order to be able to isolate shorter thrust arcs more
accurately, we set $\rho_{\text{min}}=1.0\times 10^{-6}$, which shows a small
increase in the thrust magnitude at which this bifurcation occurs. The blown-
up part shows the time when the switching function (blue line) nearly touches
the $S=0$ line, whereas the red line shows the derivative of the switching
function with respect to time, $\dot{S}$. Figure 78 shows the trajectory where
the point of the creation of the new thrust arc is shown. The final phase of
the flight consists of a thrust arc where the trajectory undergoes a sudden
turn an reaches the entry point to the L1 orbit with a final mass of
$m_{f}=1394.85$ kg. The trajectory looks different when plotted in a non-
rotating frame as is shown in Figure 79, where the blue line denotes the path
made by the entry point along the time of flight. It also shows a 3D view of
the trajectory in which the out-of-plane motion is noticeable.

Figure 79: Inertial trajectory for GTO-to-Halo problem with $T\approx 24.65$
N.

##  8 Computational Effort

In this section, computational effort is discussed. The proposed methodology
relies on solving TPBVPs to generate switching functions, which, as thrust is
swept, gives rise to neighboring extremals that underly the optimal switching
surface. The first step has to do with the determination of the minimum-thrust
solution for each $N_{\text{rev}}$. The second step, deals with generating the
corresponding switching surfaces. A third step, if one desires to determine
the idealized impulsive trajectories, is to perform $N$-impulse optimization.

Solving TPBVPs efficiently and reliably in the presence of high nonlinearity,
and high dimensionality is a challenging task [109, 110] and there are several
approaches to enhance convergence when discontinuous inputs are forcing the
differential equations [71]. Our experience is that for a great variety of
optimal trajectory cases, the hyperbolic tangent smoothing is quite helpful
and alleviates a number of difficulties associated with solving the TPBVPs
associated with switching-type controls [111]. On the other hand, if our goal
is to determine a set of impulsive solutions, in order to generate the results
similar to those reported in Table 11, it is not (usually) required to sweep
over the thrust magnitudes, in a dense incremental fashion. Instead, it is
possible to choose a relatively large thrust level and solve the corresponding
minimum-fuel TPBVP directly. Even though, we use the hyperbolic tangent
smoothing, only the final trajectory needs to approach the very sharp control
switches. As mentioned earlier, the impulsive solution reveals itself at high-
thrust levels encountered in the family of optimal trajectory problems under
consideration.

This strategy has been adopted to generate the results in Table 11. For the
Earth-to-Dionysus problem, the computation wall-time to obtain an impulsive
solution is usually only 30 to 60 seconds. This time takes into account the
time it takes to solve minimum-thrust problem according to the methodology
outlined in in this paper while relying on hyperbolic tangent smoothing,
converging with a large value of thrust that solves the corresponding minimum-
fuel TPBVP, including the optimal switching function, and the time it takes to
solve the $N$-impulse optimization problem. The $N$-impulse optimization is
the least computationally demanding part of the overall scheme. Of course,
harder problems with highly nonlinear trajectories and more impulses require
more computation time. The GTO-to-GEO orbit transfer problem proved to be the
most computationally demanding case when two-body dynamics was considered.
Generating the GTO-to-GEO impulsive solutions took several minutes.

All computations were performed on a Core i7 desktop machine with two 3.4-GHz
processors and 16 GB of RAM and the propagation of the differential equations
was achieved using a compiled version of MATLAB’s ode45 function to speed up
the numerical simulations. The absolute and relative integration convergence
tolerances were set to $1.0\times 10^{-10}$ in all simulations. For MATLAB’s
fmincon optimizer, the function and constraint tolerances were set to
$TolFun=10\times 10^{-6}$ and $TolCon=10\times 10^{-7}$, respectively.

##  9 Conclusion

We used indirect methods of optimal control theory along with Lawden’s primer
vector theory to introduce illuminating switching surfaces for large and low-
thrust minimum-fuel trajectories. These surfaces establish the connection
between impulsive and continuous-thrust trajectories. The impulsive limits
tell us “how many impulses” are associated with each application of
$N_{\text{rev}}$, the number of revolutions.

We have shown that there exists a fundamental minimum-thrust solution that
plays a pivotal role in determining fundamental switching surface, and an
associated $N^{*}_{\text{rev}}$, which in turn reveals the “optimal” number of
thrust arcs for any finite specification of maximum thrust. The fundamental
switching surface is generated by a homotopic sweep of the maximum thrust away
from the the minimum-thrust extremal among all minimum-thrust solutions,
considering all feasible multi-revolution solutions. This fundamental minimum-
thrust solution is used to construct the switching surfaces that reveal the
number and approximate time, direction and magnitudes of the impulses
associated with the optimal impulsive solution as a by-product. A number of
optimal orbit transfers with different number of revolutions for the
fundamental extremal are studied and the corresponding impulsive solution with
as many as eleven impulses are found, which demonstrate the capability of the
proposed construct. According to these results and our experience, the optimal
number of thrusting arcs at high thrust levels depend strongly on the size and
shape of the initial and final orbits, time of flight, angular phasing and the
orbits’ relative orientation.

While the switching surface can be generated for very high thrust values to
approach with arbitrary precision the impulsive limit, we find that after the
elapsed thrust-on time of all short thrust arcs is less than some tolerance
(e.g., $<0.001\times$ (time of flight)), then the impulsive approximation can
be invoked with a good accuracy. A direct method can then be initiated
accurately to converge to the impulsive solution and satisfy the force model
and all boundary conditions to a high precision. It is vital to note that when
mass variation is accounted for, considering any real propulsion system, that
minimizing $\Delta v$ does not minimize fuel consumption (or maximize payload
mass for a given engine). We showed two examples where minimizing $\Delta v$
led to multiple distinct local extremals, with identical minimum $\Delta v$
values; the local extremals $\Delta v$ matched to 7 digits! Only one of these
extremals belonged to the fundamental switching surface that minimizes fuel
consumption. So, we conclude that Edelbaum’s question does not have a unique
answer, in general. However, invoking the minimum propellant consumption for a
specific engine model gives a unique extremal.

We showed that other secondary criteria, such as time of flight (for cases
with early arrival) can also be invoked, for the same total $\Delta v$ and
modest fuel increases. The results clearly indicate the intuitively reasonable
truth that the optimal short thrusting arcs as well as limiting case of
impulses are predominantly applied near peripasis/apoapsis and/or the
ascending or descending nodes, as might be anticipated for orbit
raising/lowering and changes in inclination.

Finally, we note that the switching surfaces provide a unified means for
optimizing low thrust, high thrust, and impulsive maneuvers, and enable
important global insights in mission design. In principle, the methodology
proposed in this paper can be pursued for many engineering systems if/when a
systematic study of the resulting switching surfaces are performed by sweeping
various important parameters of interest.

## Conflict of Interest Statement

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## Acknowledgments

We are pleased like to thank our sponsors: AFOSR (Stacie Williams), AFRL (Alok
Das et al.) and ADS (Matt Wilkins), for their support and collaborations under
various contracts and grants. We are indebted to many colleagues, we
especially appreciate the interactions we have had on the ideas in this paper
with Manoranjan Majji.

## References

  * [1] Petropoulos, A. E. and Longuski, J. M., “Shape-based algorithm for the automated design of low-thrust, gravity assist trajectories,” Journal of Spacecraft and Rockets, Vol. 41, No. 5, 2004, pp. 787–796, doi:10.2514/1.13095.
  * [2] Taheri, E. and Abdelkhalik, O., “Shape based approximation of constrained low-thrust space trajectories using Fourier series,” Journal of Spacecraft and Rockets, Vol. 49, No. 3, 2012, pp. 535–546, doi:10.2514/1.58789.
  * [3] Zuiani, F., Vasile, M., Palmas, A., and Avanzini, G., “Direct transcription of low-thrust trajectories with finite trajectory elements,” Acta Astronautica, Vol. 72, 2012, pp. 108–120, doi:10.1016/j.actaastro.2011.09.011.
  * [4] Avanzini, G., Palmas, A., and Vellutini, E., “Solution of low-thrust lambert problem with perturbative expansions of equinoctial elements,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 9, 2015, pp. 1585–1601, doi:10.2514/1.G001018.
  * [5] Kluever, C. A., “Designing Transfers to Geostationary Orbit Using Combined Chemical–Electric Propulsion,” Journal of Spacecraft and Rockets, Vol. 52, No. 4, 2015, pp. 1144–1151, doi:10.2514/1.A33259.
  * [6] Taheri, E. and Abdelkhalik, O., “Initial three-dimensional low-thrust trajectory design,” Advances in Space Research, Vol. 57, No. 3, 2016, pp. 889–903, doi:10.1016/j.asr.2015.11.034.
  * [7] Zeng, K., Geng, Y., and Wu, B., “Shape-based analytic safe trajectory design for spacecraft equipped with low-thrust engines,” Aerospace Science and Technology, Vol. 62, 2017, pp. 87–97, doi:10.1016/j.ast.2016.12.006.
  * [8] Roa, J., Peláez, J., and Senent, J., “New analytic solution with continuous thrust: generalized logarithmic spirals,” Journal of Guidance, Control, and Dynamics, Vol. 39, No. 10, 2016, pp. 2336–2351, doi:10.2514/1.G000341.
  * [9] Taheri, E., Kolmanovsky, I., and Atkins, E., “Shaping low-thrust trajectories with thrust-handling feature,” Advances in Space Research, Vol. 61, 2018, pp. 879–890, doi:10.1016/j.asr.2017.11.006.
  * [10] Longuski, J. M. and Williams, S. N., “Automated design of gravity-assist trajectories to Mars and the outer planets,” Celestial Mechanics and Dynamical Astronomy, Vol. 52, No. 3, 1991, pp. 207–220, doi:10.1007/BF00048484.
  * [11] Izzo, D., Becerra, V. M., Myatt, D. R., Nasuto, S. J., and Bishop, J. M., “Search space pruning and global optimisation of multiple gravity assist spacecraft trajectories,” Journal of Global Optimization, Vol. 38, No. 2, 2007, pp. 283–296, doi:10.1007/s10898–006–9106–0.
  * [12] Olympio, J. T. and Izzo, D., “Designing optimal multi-gravity-assist trajectories with free number of impulses,” International Symposium on Space Flights Dynamics, Toulouse, France. ESA ESTEC, 2009.
  * [13] Vasile, M. and Ceriotti, M., “8 Incremental Techniques for Global Space Trajectory Design,” Spacecraft Trajectory Optimization, edited by B. Conway, Vol. 29, 2010, pp. 202–237, doi:10.1017/CBO9780511778025.
  * [14] Englander, J. A., Conway, B. A., and Williams, T., “Automated mission planning via evolutionary algorithms,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 6, 2012, pp. 1878–1887, doi:10.2514/1.54101.
  * [15] Abdelkhalik, O. and Gad, A., “Dynamic-size multiple populations genetic algorithm for multigravity-assist trajectory optimization,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 2, 2012, pp. 520–529, doi:10.2514/1.54330.
  * [16] Chilan, C. M. and Conway, B. A., “Automated design of multiphase space missions using hybrid optimal control,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 5, 2013, pp. 1410–1424, doi:.
  * [17] Ellison, D. H., Conway, B. A., Englander, J. A., and Ozimek, M. T., “Analytic Gradient Computation for Bounded-Impulse Trajectory Models Using Two-Sided Shooting,” Journal of Guidance, Control, and Dynamics, Vol. 41, No. 7, 2018, pp. 1449–1462, doi:10.2514/1.G003077.
  * [18] Landau, D., “Efficient Maneuver Placement for Automated Trajectory Design,” Journal of Guidance, Control, and Dynamics, Vol. 41, No. 7, 2018, pp. 1531–1541, doi:10.2514/1.G003172.
  * [19] Edelbaum, T. N., “How many impulses?” Astronautics and Aeronautics, Vol. 5, No. 11, 1967, pp. 64–69, doi:10.2514/6.1966–7.
  * [20] Sims, J. A. and Flanagan, S. N., “Preliminary design of low-thrust interplanetary missions,” AAS/AIAA Astrodynamics Specialist Conference, AAS Paper 99-338, Girdwood, Alaska,, Aug. 1997.
  * [21] Gooding, R., “A procedure for the solution of Lambert’s orbital boundary-value problem,” Celestial Mechanics and Dynamical Astronomy, Vol. 48, No. 2, 1990, pp. 145–165.
  * [22] Battin, R. H., An introduction to the mathematics and methods of astrodynamics, AIAA, 1999, chp. 7.
  * [23] Izzo, D., “Revisiting Lambert’s problem,” Celestial Mechanics and Dynamical Astronomy, Vol. 121, No. 1, 2015, pp. 1–15, doi:10.1007/s10569–014–9587–y.
  * [24] Arora, N., Russell, R. P., Strange, N., and Ottesen, D., “Partial derivatives of the solution to the Lambert boundary value problem,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 9, 2015, pp. 1563–1572, doi:10.2514/1.G001030.
  * [25] Bombardelli, C., Gonzalo, J. L., and Roa, J., “Approximate analytical solution of the multiple revolution Lambert’s targeting problem,” Journal of Guidance, Control, and Dynamics, Vol. 41, No. 3, 2018, pp. 792–801, doi:10.2514/1.G002887.
  * [26] Vaddi, S., Alfriend, K. T., Vadali, S., and Sengupta, P., “Formation establishment and reconfiguration using impulsive control,” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 2, 2005, pp. 262–268, doi:10.2514/1.6687.
  * [27] Gurfil, P., “Relative motion between elliptic orbits: generalized boundedness conditions and optimal formation keeping,” Journal of Guidance, Control, and Dynamics, Vol. 28, No. 4, 2005, pp. 761–767, doi:10.2514/1.9439.
  * [28] Woffinden, D. C. and Geller, D. K., “Navigating the road to autonomous orbital rendezvous,” Journal of Spacecraft and Rockets, Vol. 44, No. 4, 2007, pp. 898–909, doi:10.2514/1.30734.
  * [29] Ichikawa, A. and Ichimura, Y., “Optimal impulsive relative orbit transfer along a circular orbit,” Journal of guidance, control, and dynamics, Vol. 31, No. 4, 2008, pp. 1014–1027, doi:10.2514/1.32820.
  * [30] Kim, D.-Y., Woo, B., Park, S.-Y., and Choi, K.-H., “Hybrid optimization for multiple-impulse reconfiguration trajectories of satellite formation flying,” Advances in Space Research, Vol. 44, No. 11, 2009, pp. 1257–1269, doi:10.1016/j.asr.2009.07.029.
  * [31] Huang, W., “Optimal multi-impulse orbit transfer using nonlinear relative motion dynamics,” The Journal of the Astronautical Sciences, Vol. 59, No. 1-2, 2012, pp. 237–258, doi:10.1007/s40295–013–0015–4.
  * [32] Spencer, D. A. and Lovell, T. A., “Maneuver Design Using Relative Orbital Elements,” The Journal of the Astronautical Sciences, Vol. 62, No. 4, 2015, pp. 315–350, doi:10.1007/s40295–015–0072–y.
  * [33] Roscoe, C. W., Westphal, J. J., Griesbach, J. D., and Schaub, H., “Formation Establishment and Reconfiguration Using Differential Elements in J 2-Perturbed Orbits,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 9, 2015, pp. 1725–1740, doi:10.2514/1.G000999.
  * [34] Sobiesiak, L. A. and Damaren, C. J., “Impulsive Spacecraft Formation Maneuvers with Optimal Firing Times,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 10, 2015, pp. 1994–2000, doi:10.2514/1.G001095.
  * [35] Starek, J. A., Schmerling, E., Maher, G. D., Barbee, B. W., and Pavone, M., “Fast, Safe, Propellant-Efficient Spacecraft Motion Planning Under Clohessy–Wiltshire–Hill Dynamics,” Journal of Guidance, Control, and Dynamics, Vol. 40, No. 2, 2016, pp. 418–438, doi:10.2514/1.G001913.
  * [36] Chernick, M. and D’Amico, S., “New closed-form solutions for optimal impulsive control of spacecraft relative motion,” Journal of Guidance, Control, and Dynamics, Vol. 41, No. 2, 2018, pp. 301–319, doi:10.2514/1.G002848.
  * [37] Vinh, N. X., Gilbert, E. G., Howe, R. M., Sheu, D., and Lu, P., “Reachable domain for interception at hyperbolic speeds,” Acta Astronautica, Vol. 35, No. 1, 1995, pp. 1–8, doi:10.1016/0094–5765(94)00132–6.
  * [38] Xue, D., Li, J., Baoyin, H., and Jiang, F., “Reachable domain for spacecraft with a single impulse,” Journal of guidance, control, and dynamics, Vol. 33, No. 3, 2010, pp. 934–942, doi:10.2514/1.43963.
  * [39] Zhang, G., Cao, X., and Ma, G., “Reachable domain of spacecraft with a single tangent impulse considering trajectory safety,” Acta Astronautica, Vol. 91, 2013, pp. 228–236, doi:10.1016/j.actaastro.2013.06.016.
  * [40] Wen, C., Zhao, Y., Shi, P., and Hao, Z., “Orbital accessibility problem for spacecraft with a single impulse,” Journal of Guidance, Control, and Dynamics, Vol. 37, No. 4, 2014, pp. 1260–1271, doi:10.2514/1.62629.
  * [41] Holzinger, M. J., Scheeres, D. J., and Erwin, R. S., “On-Orbit Operational Range Computation Using Gauss’s Variational Equations with J 2 Perturbations,” Journal of Guidance, Control, and Dynamics, Vol. 37, No. 2, 2014, pp. 608–622, doi:10.2514/1.53861.
  * [42] Wen, C., Zhao, Y., and Shi, P., “Precise determination of reachable domain for spacecraft with single impulse,” Journal of Guidance, Control, and Dynamics, Vol. 37, No. 6, 2014, pp. 1767–1779, doi:10.2514/1.G000583.
  * [43] Chen, Q., Qiao, D., Shang, H., and Liu, X., “A new method for solving reachable domain of spacecraft with a single impulse,” Acta Astronautica, Vol. 145, 2018, pp. 153–164, doi:10.1016/j.actaastro.2018.01.040.
  * [44] Lawden, D. F., Optimal trajectories for space navigation, Butterworths, 1963, chp.3.
  * [45] Handelsman, M. and Lion, P., “Primer vector on fixed-time impulsive trajectories.” AIAA Journal, Vol. 6, No. 1, 1968, pp. 127–132, doi:10.2514/3.4452.
  * [46] Jezewski, D. J. and Rozendaal, H. L., “An efficient method for calculating optimal free-space n-impulse trajectories.” AIAA journal, Vol. 6, No. 11, 1968, pp. 2160–2165, doi:10.2514/3.4949.
  * [47] Prussing, J. E., “2 Primer Vector Theory and Applications,” Spacecraft Trajectory Optimization, edited by B. Conway, Vol. 29, 2010, pp. 16–36, doi:10.1017/CBO9780511778025.
  * [48] Prussing, J. E. and Conway, B. A., Orbital mechanics, Oxford University Press, USA, 1993.
  * [49] Prussing, J. E., “A class of optimal two-impulse rendezvous using multiple-revolution Lambert solutions,” The Journal of the Astronautical Sciences, Vol. 48, No. 2-3, 2000, pp. 131–148.
  * [50] Prussing, J. E. and Chiu, J.-H., “Optimal multiple-impulse time-fixed rendezvous between circular orbits,” Journal of Guidance, Control, and Dynamics, Vol. 9, No. 1, 1986, pp. 17–22, doi:10.2514/3.20060.
  * [51] Hughes, S. P., Mailhe, L. M., and Guzman, J. J., “A comparison of trajectory optimization methods for the impulsive minimum fuel rendezvous problem,” 2003.
  * [52] Arzelier, D., Bréhard, F., Deak, N., Joldes, M., Louembet, C., Rondepierre, A., and Serra, R., “Linearized impulsive fixed-time fuel-optimal space rendezvous: A new numerical approach,” IFAC-PapersOnLine, Vol. 49, No. 17, 2016, pp. 373–378, doi:10.1016/j.ifacol.2016.09.064.
  * [53] Abdelkhalik, O. and Mortari, D., “N-impulse orbit transfer using genetic algorithms,” Journal of Spacecraft and Rockets, Vol. 44, No. 2, 2007, pp. 456–460, doi:10.2514/1.24701.
  * [54] Pontani, M., Ghosh, P., and Conway, B. A., “Particle swarm optimization of multiple-burn rendezvous trajectories,” Journal of Guidance, Control, and Dynamics, Vol. 35, No. 4, 2012, pp. 1192–1207, doi:10.2514/1.55592.
  * [55] Luo, Y.-Z., Tang, G.-J., Li, Y.-J., and Li, H.-Y., “Optimization of multiple-impulse, multiple-revolution, rendezvous-phasing maneuvers,” Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4, 2007, pp. 946–952, doi:10.2514/1.25620.
  * [56] Sentinella, M. R. and Casalino, L., “Cooperative evolutionary algorithm for space trajectory optimization,” Celestial Mechanics and Dynamical Astronomy, Vol. 105, No. 1-3, 2009, pp. 211.
  * [57] Luo, Y.-Z., Zhang, J., Li, H.-y., and Tang, G.-J., “Interactive optimization approach for optimal impulsive rendezvous using primer vector and evolutionary algorithms,” Acta Astronautica, Vol. 67, No. 3-4, 2010, pp. 396–405, doi:10.1016/j.actaastro.2010.02.014.
  * [58] Colasurdo, G. and Pastrone, D., “Indirect optimization method for impulsive transfers,” Astrodynamics Conference, 1994, p. 3762.
  * [59] Shen, H.-X., Casalino, L., and Luo, Y.-Z., “Global search capabilities of indirect methods for impulsive transfers,” The Journal of the Astronautical Sciences, Vol. 62, No. 3, 2015, pp. 212–232, doi:10.1007/s40295–015–0073–x.
  * [60] Neustadt, L. W., “Optimization, a moment problem, and nonlinear programming,” Journal of the Society for Industrial and Applied Mathematics, Series A: Control, Vol. 2, No. 1, 1964, pp. 33–53, doi:10.1137/0302004.
  * [61] Prussing, J. E., “Optimal impulsive linear systems: Sufficient conditions and maximum number of impulses,” The Journal of the Astronautical Sciences, Vol. 43, No. 2, 1995, pp. 195–206, doi:.
  * [62] Carter, T., “Optimal impulsive space trajectories based on linear equations,” Journal of optimization theory and applications, Vol. 70, No. 2, 1991, pp. 277–297, doi:10.1007/BF00940627.
  * [63] Carter, T. and Brient, J., “Fuel-optimal rendezvous for linearized equations of motion,” Journal of Guidance, Control, and Dynamics, Vol. 15, No. 6, 1992, pp. 1411–1416, doi:.
  * [64] Lawden, D. F., “Optimal transfers between coplanar elliptical orbits,” Journal of guidance, control, and dynamics, Vol. 15, No. 3, 1992, pp. 788–791, doi:10.2514/3.20909.
  * [65] Baranov, A. and Roldugin, D., “Six-impulse maneuvers for rendezvous of spacecraft in near-circular noncoplanar orbits,” Cosmic Research, Vol. 50, No. 6, 2012, pp. 441–448, doi:10.1134/S0010952512050012.
  * [66] Iorfida, E., Palmer, P., and Roberts, M., “Geometric approach to the perpendicular thrust case for trajectory optimization,” Journal of Guidance, Control, and Dynamics, 2016, pp. 1059–1068, doi:10.2514/1.G001525.
  * [67] Serra, R., Arzelier, D., and Rondepierre, A., “Analytical solutions for impulsive elliptic out-of-plane rendezvous problem via primer vector theory,” IEEE Transactions on Control Systems Technology, Vol. 26, No. 1, 2018, pp. 207–221, doi:10.1109/TCST.2017.2656022.
  * [68] Neustadt, L. W., “A general theory of minimum-fuel space trajectories,” Journal of the Society for Industrial and Applied Mathematics, Series A: Control, Vol. 3, No. 2, 1965, pp. 317–356, doi:10.1137/0303023.
  * [69] Gergaud, J. and Haberkorn, T., “Orbital transfer: some links between the low-thrust and the impulse cases,” Acta Astronautica, Vol. 60, No. 8-9, 2007, pp. 649–657, doi:10.1016/j.actaastro.2006.10.009.
  * [70] Zhu, Z., Gan, Q., Yang, X., and Gao, Y., “Solving fuel-optimal low-thrust orbital transfers with bang-bang control using a novel continuation technique,” Acta Astronautica, Vol. 137, 2017, pp. 98–113, doi:10.1016/j.actaastro.2017.03.032.
  * [71] Taheri, E., Kolmanovsky, I., and Atkins, E., “Enhanced smoothing technique for indirect optimization of minimum-fuel low-thrust trajectories,” Journal of Guidance, Control, and Dynamics, Vol. 39, No. 11, 2016, pp. 2500–2511, doi:10.2514/1.G000379.
  * [72] Kluever, C. A. and Pierson, B. L., “Optimal earth-moon trajectories using nuclear electric propulsion,” Journal of Guidance Control and Dynamics, Vol. 20, No. 2, 1997, pp. 239–245, doi:10.2514/2.4058.
  * [73] Caillau, J.-B. and Noailles, J., “Coplanar control of a satellite around the Earth,” ESAIM: Control, Optimisation and Calculus of Variations, Vol. 6, 2001, pp. 239–258.
  * [74] Bertrand, R. and Epenoy, R., “New smoothing techniques for solving bang–bang optimal control problems—numerical results and statistical interpretation,” Optimal Control Applications and Methods, Vol. 23, No. 4, 2002, pp. 171–197, doi:10.1002/oca.709.
  * [75] Caillau, J.-B., Gergaud, J., and Noailles, J., “3D geosynchronous transfer of a satellite: Continuation on the thrust,” Journal of optimization theory and applications, Vol. 118, No. 3, 2003, pp. 541–565, doi:10.1023/B:JOTA.0000004870.74778.ae.
  * [76] Haberkorn, T., Martinon, P., Gergaud, J., et al., “Low-thrust minimum-fuel orbital transfer: a homotopic approach,” Journal of Guidance Control and Dynamics, Vol. 27, 2004, pp. 1046–1060, doi:10.2514/1.4022.
  * [77] Bonnard, B., Caillau, J.-B., and Trélat, E., “Geometric optimal control of elliptic Keplerian orbits,” Discrete and Continuous Dynamical Systems-Series S, 2005, pp. 929–956.
  * [78] La Mantia, M. and Casalino, L., “Indirect optimization of low-thrust capture trajectories,” Journal of Guidance Control and Dynamics, Vol. 29, No. 4, 2006, pp. 1011–1014, doi:10.2514/1.18986.
  * [79] Russell, R. P., “Primer vector theory applied to global low-thrust trade studies,” Journal of Guidance Control and Dynamics, Vol. 30, No. 2, 2007, pp. 460–472, doi:10.2514/1.22984.
  * [80] Silva, C. and Trélat, E., “Smooth regularization of bang-bang optimal control problems,” IEEE Transactions on Automatic Control, Vol. 55, No. 11, 2010, pp. 2488–2499, doi:10.1109/TAC.2010.2047742.
  * [81] Jiang, F., Baoyin, H., and Li, J., “Practical techniques for low-thrust trajectory optimization with homotopic approach,” Journal of Guidance, Control and Dynamics, Vol. 35, No. 1, 2012, pp. 245–258, doi:10.2514/1.52476.
  * [82] Dutta, A., Arora, N., and Russell, R. P., “Peer-to-peer refueling strategy using low-thrust propulsion,” Journal of Spacecraft and Rockets, Vol. 49, No. 5, 2012, pp. 944–954, doi:10.2514/1.A32106.
  * [83] Pan, B., Lu, P., Pan, X., and Ma, Y., “Double-homotopy method for solving optimal control problems,” Journal of Guidance, Control, and Dynamics, Vol. 39, No. 8, 2016, pp. 1706–1720,doi:10.2514/1.G001553.
  * [84] Zhao, S. and Zhang, J., “Minimum-fuel station-change for geostationary satellites using low-thrust considering perturbations,” Acta Astronautica, Vol. 127, 2016, pp. 296–307, doi:10.1016/j.actaastro.2016.05.028.
  * [85] Zhao, S., Gurfil, P., and Zhang, J., “Initial Costates for Low-Thrust Minimum-Time Station Change of Geostationary Satellites,” Journal of Guidance, Control, and Dynamics, Vol. 39, No. 12, 2016, pp. 2746–2756, doi:10.2514/1.G000431.
  * [86] Chi, Z., Yang, H., Chen, S., and Li, J., “Homotopy method for optimization of variable-specific-impulse low-thrust trajectories,” Astrophysics and Space Science, Vol. 362:216, 2017, doi:10.1007/s10509-017-3196-7.
  * [87] Taheri, E., Li, N. I., and Kolmanovsky, I., “Co-state initialization for the minimum-time low-thrust trajectory optimization,” Advances in Space Research, Vol. 59, No. 9, 2017, pp. 2360–2373, doi:10.1016/j.asr.2017.02.010.
  * [88] Mall, K. and Grant, M. J., “Epsilon-Trig Regularization Method for Bang-Bang Optimal Control Problems,” Journal of Optimization Theory and Applications, Vol. 174, No. 2, 2017, pp. 500–517, doi:10.1007/s10957–017–1129–9.
  * [89] Sullo, N., Peloni, A., and Ceriotti, M., “Low-thrust to solar-sail trajectories: a homotopic approach,” Journal of Guidance, Control, and Dynamics, Vol. 40, No. 11, 2017, pp. 2796–2806, doi:10.2514/1.G002552.
  * [90] Zhu, J., Trélat, E., and Cerf, M., “Geometric optimal control and applications to aerospace,” Pacific Journal of Mathematics for Industry, Vol. 9, No. 1, 2017, pp. 8, doi:10.1186/s40736–017–0033–4.
  * [91] Saranathan, H. and Grant, M. J., “Relaxed Autonomously Switched Hybrid System Approach to Indirect Multiphase Aerospace Trajectory Optimization,” Journal of Spacecraft and Rockets, Vol. 55, No. 3, 2017, pp. 611–621, doi:10.2514/1.A34012.
  * [92] Saghamanesh, M. and Baoyin, H., “A robust homotopic approach for continuous variable low-thrust trajectory optimization,” Advances in Space Research, Vol. 62, No. 11, 2018, pp. 3095–3113, doi:10.1016/j.asr.2018.08.046.
  * [93] Pérez-Palau, D. and Epenoy, R., “Fuel optimization for low-thrust Earth–Moon transfer via indirect optimal control,” Celestial Mechanics and Dynamical Astronomy, Vol. 130, No. 2, 2018, pp. 21, doi:10.1007/s105.
  * [94] Junkins, J. L. and Taheri, E., “Exploration of Alternative State Vector Choices for Low-Thrust Trajectory Optimization,” Journal of Guidance, Control, and Dynamics, Vol. 42, No. 1, 2018, pp. 47–64, doi:10.2514/1.G003686.
  * [95] Pan, B., Pan, X., and Zhang, S., “A new probability-one homotopy method for solving minimum-time low-thrust orbital transfer problems,” Astrophysics and Space Science, Vol. 363, No. 9, 2018, pp. 198, doi:10.1007/s10509–018–3420–0.
  * [96] Taheri, E. and Junkins, J. L., “Generic Smoothing for Optimal Bang-Off-Bang Spacecraft Maneuvers,” Journal of Guidance, Control, and Dynamics, Vol. 41, No. 11, 2018, pp. 2470–2475, doi:10.2514/1.G003604.
  * [97] Pan, B., Pan, X., and Ma, Y., “A quadratic homotopy method for fuel-optimal low-thrust trajectory design,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, Vol. 233, No. 5, 2019, pp. 1741–1757, doi:10.1177/0954410018761965.
  * [98] Walker, M., “A set of modified equinoctial orbit elements,” Celestial Mechanics and Dynamical Astronomy, Vol. 38, No. 4, 1986, pp. 391–392, doi:10.1007/BF01227493.
  * [99] Taheri, E. and Junkins, J., “A Generic Approach for Optimal Bang-Off-Bang Spacecraft Maneuvers,” 41st Annual AAS Guidance & Control Conference, Breckenridge, Colorado, 2018, AAS 18-088.
  * [100] Kelley, H. J., “A second variation test for singular extremals,” AIAA Journal, Vol. 2, No. 8, 1964, pp. 1380–1382, doi:10.2514/3.2562.
  * [101] Breakwell, J., “Minimum impulse transfer,” Astrodynamics Conference, 1964, p. 416.
  * [102] Small, H. W., “Minimum N-impulse time-free transfer between elliptic orbits,” AIAA Journal, Vol. 9, No. 4, 1971, pp. 594–599.
  * [103] Small, H. W., “Globally optimal parking orbit transfer,” Astrodynamics Conference, 1982, p. 1478.
  * [104] Prussing, J., “Optimal two-and three-impulse fixed-time rendezvous in the vicinity of a circular orbit,” AIAA Journal, Vol. 8, No. 7, 1970, pp. 1221–1228, doi:10.2514/2.7041.
  * [105] Woollands, R. M., Read, J. L., Probe, A. B., and Junkins, J. L., “Multiple Revolution Solutions for the Perturbed Lambert Problem using the Method of Particular Solutions and Picard Iteration,” The Journal of the Astronautical Sciences, Vol. 64, No. 4, 2017, pp. 361–378, doi:10.1007/s40295–017–0116–6.
  * [106] Taheri, E. and Abdelkhalik, O., “Fast initial trajectory design for low-thrust restricted-three-body problems,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 11, 2015, pp. 2146–2160, doi:10.2514/1.G000878.
  * [107] Edelbaum, T. N., “Propulsion requirements for controllable satellites,” ARS Journal, Vol. 31, No. 8, 1961, pp. 1079–1089, doi:10.2514/8.5723.
  * [108] Vinh, N. X., Kuo, S. H., and Marchal, C., “Optimal time-free nodal transfers between elliptical orbits,” Acta Astronautica, Vol. 17, No. 8, 1988, pp. 875–880, doi:10.1016/0094–5765(88)90170–1.
  * [109] Betts, J. T., “Survey of numerical methods for trajectory optimization,” Journal of Guidance Control and Dynamics, Vol. 21, No. 2, 1998, pp. 193–207, doi:10.2514/2.4231.
  * [110] Woolley, R. C., Baker, J. D., Landau, D. F., and Nicholas, A. K., “Cargo logistics for a Notional Mars Base using Solar Electric Propulsion,” Acta Astronautica, Vol. 156, 2019, pp. 51–57, doi:10.1016/j.actaastro.2018.08.026.
  * [111] Woollands, R. M., Taheri, E., and Junkins, J. L., “Efficient Computation of Optimal Low Thrust Gravity Perturbed Orbit Transfers,” The Journal of the Astronautical Sciences, Vol. -, No. -, 2019, pp. 1–27, doi:10.1007/s40295–019–00152–9.


# Supplementary material for “Long-range influence of a pump on a critical
fluid”

##  I Numerical methods

In this section we describe the numerical methods of this work.

Kawasaki Dynamics: To simulate directly the conserved dynamics described in
the main text, we use standard local Kawasaki dynamics on a two-dimensional
Ising lattice gas of dimensions $L\times L$. The system is initiated with a
prescribed total magnetization by setting the relative number of $s=+1$ and
$s=-1$ spins. At each Monte Carlo step, a pair of nearest-neighbor spins is
chosen at random. If the spins are opposite, they are exchanged with
probability $p_{\rm ex}=\min\left[\exp\left(-\beta\Delta E\right),1\right]$
where $\beta$ is the inverse temperature and $\Delta E$ is the energy
difference between the exchanged and initial states which is obtained from the
Ising energy $E(\\{s_{i}\\})=-J\sum_{\langle i,j\rangle}s_{i}s_{j}$. Here
$J>0$ is the coupling constant and $\langle i,j\rangle$ denotes nearest
neighbor sites. To enhance the current and density modulation the pumping acts
along two bonds at the center of the lattice – a bond between sites
$(L/2,L/2)$ and $(L/2+1,L/2)$ and another bond between sites $(L/2,L/2+1)$ and
$(L/2+1,L/2+1)$. The updating rule of the two pump bonds is different. If
$s_{i}=+1$ on the left site and $s_{i}=-1$ on the right site of either pair of
sites they are exchanged. This corresponds to an infinite biasing rate. The
possible transitions are depicted in Fig. 1.

The Kawasaki dynamics was used to study the density profile below the critical
point and the currents in the system across the temperature range, described
in the SM below. The current along each bond is measured by counting the
exchange events per unit time and the streamlines appearing in Fig. 3 of the
SM were obtained using a standard Matlab tool. In addition, the Kawasaki
algorithm was used to test the agreement between the steady-state density with
that predicted by the effective equilibrium model.

In panels (j)-(l) of Figure 1 of the main text, the density profile obtained
by numerical simulation below the critical point is displayed. It was carried
out on an Ising square lattice of size $L\times L$ with $L=256$, imposing
helical boundary conditions [1]. This means that the last spin in each row is
coupled to the first spin of the next row instead of the first spin of the
same row. This is computationally convenient as it allows one to store the
lattice in a linear array. To do so we number the lattice sites from $1$ to
$L^{2}$ along the rows, with the first row labelled $1$ to $L$. Then the
adjacency rule has a simple description in which the neighbors of the $i^{\rm
th}$ spin on the two-dimensional lattice are situated at $(i\pm 1)$ modulo
$L^{2}$ and at $(i\pm L)$ modulo $L^{2}$. The density profile was averaged
over 300 realizations, running for $1.5\cdot 10^{6}$ Monte Carlo sweeps.

Figure 1: An illustration of the possible exchanges on the two-dimensional
lattice, here represented as a lattice gas with $s_{i}=+1$ for an occupied
site and $s_{i}=-1$ for an empty site. All arrows denote possible exchanges of
an occupied site moving into an empty site. Blue arrows stand for transitions
associated with the pump. In our numerics these occur whenever an exchange is
possible.

Wolff’s algorithm: Away from the critical point, the system is easily
simulated using Kawasaki dynamics. However, most of our work is dedicated to
the critical behavior of the system, where it is subject to critical slowing
down [2]. Indeed, attempts to simulate large systems using Kawasaki dynamics
proved to converge extremely slowly, as expected.

To obtain the density profiles in reasonable time, we use the fact that as
long as $T$ is not far below the critical temperature, the density profile can
be obtained from a mapping onto an effective equilibrium system, as discussed
in the main text. The pump used in the Kawasaki dynamics is then replaced by a
long-ranged effective magnetic field of strength $f$. This furthermore allows
us, through equivalence of ensembles, to use non-conserving dynamics to reach
the same equilibrium density profile as that obtained using conserving
dynamics. Therefore, equilibrium, non-conserving, cluster algorithms can be
used to measure the desired density profiles at criticality. In our case, we
used an adaptation of Wolff’s algorithm [3] that allows for a (spatially-
varying) magnetic field. This is done by using a “ghost spin” – an additional
degree of freedom that is assigned to the magnetic field itself, as detailed
in Ref. [4] and outlined below.

The energy function used in the simulations is given by

$E(\\{s_{i}\\},s_{g})=-J\sum_{\langle
i,j\rangle}s_{i}s_{j}-\sum_{i}H_{i}s_{g}s_{i}\;,$ (1)

with $s_{i},s_{j}\in\\{-1,1\\}$ and $H_{i}$ is the value of the field on the
$i^{\rm th}$ site. The difference with the ordinary Ising Hamiltonian is the
addition of the ghost spin $s_{g}\in\\{-1,1\\}$ which couples to all the spins
in the system and sets the global direction of the magnetic field. With it,
the $\mathbb{Z}_{2}$ symmetry of the system in the absence of a magnetic field
is restored, at the expense of measuring the state of the system relative to
the fluctuating direction set by the ghost spin. For completeness, we briefly
outline the algorithm [4] used in our simulation. In addition to the spin
configuration the algorithm keeps track of two lists of spins. The first is
called a cluster and the second a queue. The algorithm repeatedly builds
clusters, flips them, and restarts them. The queue is used to keep track of
spins that are in the process of being added to the cluster. Each cycle is
carried out as follows:

  1. 1.

Start with an empty cluster and an empty queue. Then choose a spin (ordinary
or ghost) randomly and add it to the queue.

  2. 2.

Until the queue becomes empty:

    1. (a)

Pick the oldest spin still in the queue, say $s_{i}$.

    2. (b)

If the $i^{\rm th}$ spin $s_{i}$ is not registered in the cluster,

      1. i.

Register it to the cluster.

      2. ii.

If the $i^{\rm th}$ spin is an ordinary spin,

        1. A.

Add the ghost spin $s_{g}$ to the queue with probability
$p_{H}=\max\left[1-\exp\left(-2\beta H_{i}s_{g}s_{i}\right),0\right]$.

        2. B.

Add each of its nearest neighbors to the queue with probability
$p_{s}=1-\exp\left(-2\beta J\right)$ only if the neighboring spin has the same
spin value as the $i^{\rm th}$ spin.

      3. iii.

If the $i^{\rm th}$ spin is the ghost spin $s_{g}$ add each spin $j$ of the
lattice with probability $p_{H}=\max\left[1-\exp\left(-2\beta
H_{j}s_{g}s_{j}\right),0\right]$ to the queue.

    3. (c)

Remove the spin from the queue.

  3. 3.

After the queue becomes empty flip all the spins in the cluster and return to
1.

Finally, measuring the density involves multiplying the values of all spins by
that of the ghost spin, $s_{g}$.

In the figure of the main text, the simulations above and at criticality were
done using the Wolff algorithm presented here. The Ising system was simulated
on a $512\times 512$ square lattice, imposing helical boundary conditions, as
in the Kawasaki simulations. The temperature of the simulation above the
critical point was chosen to be $T=2T_{c}$, with $T_{c}\simeq 2.27J$. The
force strength $f$ was set to $f=1$ for $T>T_{c}$, $f=0.03$ for $T=T_{c}$
probing the near-field behavior, and $f=100$ for $T=T_{c}$ probing the far-
field behavior. The density profile above $T_{c}$ and the far-field profile at
$T_{c}$ were averaged over 250 realizations, each simulated for $10^{3}$
cycles. For the near-field density profile at criticality, the average was
done over 500 realizations, each simulated for $1.5\cdot 10^{4}$ cycles. Note
that due to the ghost spin, each cycle involves all spins in the lattice and
hence analogous to a usual Monte Carlo sweep.

Figure 2 shows the critical behavior of a three-dimensional Ising system.
These simulations were done using the same Wolff algorithm on a $64\times
64\times 64$ cubic lattice. For simplicity, here we used the magnetic field of
an infinite system (as given by Eq. (8) in the main text) and do not account
for the periodic boundary conditions. The temperature was set to the
numerically-estimated critical value for the three-dimensional cubic Ising
model [5], with $\beta_{c}\simeq 0.22165J^{-1}$ the inverse critical
temperature and $J$ the Ising coupling constant. The pump force strength $f$
was set to $f=0.5$ for probing the near-field behavior and to $f=50$ for
probing the far-field behavior. The density profiles were averaged over 100
realizations, simulated for $5\cdot 10^{3}$ cycles (analogous to a Monte Carlo
sweep) for the near-field density profile and for $10^{3}$ cycles for the far-
field profile. The results are compared with the theoretical predictions for
$r<r^{*}$ where $\phi\sim{\bf f}\cdot{\bf r}/r^{(d-2+\eta)},\,\eta\simeq
0.036$ and for $r>r^{*}$ where $\phi\sim|{\bf f}\cdot{\bf
r}/r^{d}|^{1/\delta},\,\delta\simeq 4.79$. Note that here, in contrast to two
dimensions, the theory predicts a decaying density profile for $r<r^{*}$. This
is consistent with the numerical results.

Figure 2: Results for a three dimensional $L\times L\times L$ lattice gas of
zero magnetization with $L=64$, in the presence of a pump (indicated by a red
arrow) at criticality ($T\simeq 4.51J$). The top row ((a),(d)) shows a slice
of the magnetization profiles in the plane of the pump $z=0$, the middle row
((b),(e)) the measured angular dependence of the density slice at $r=0.28L$
(marked by the dotted circles in the top row) compared to the theoretical
prediction (black line), and the bottom row ((c),(f)) the radial dependence of
the magnetization along the direction of the pump, as a function the distance
from it compared to the theoretical prediction (black line). We consider two
pump strengths allowing us to verify the behavior below and above $r^{*}$ (see
text).

##  II Current measurement

As argued in the main text a localized drive generates a current whose form is
that of a field of an electric dipole, and given by Eq. (5) above, at
criticality, and slightly below it.

To verify this, we use the Kawasaki dynamics described above. The resulting
current streamlines, as well as the $r^{-2}$ decay, at different temperatures
above, at, and below the critical point are displayed in Fig. 3. Note that
since the pump’s strength $f$ depends implicitly on the density of particles
at the origin, the amplitude of the current varies as a function of the
temperature as well. Surprisingly, while the mapping to effective equilibrium
is not expected to hold far below the critical point, due to the density
dependence of the mobility, we find that the current nonetheless exhibits the
same dipolar form as if the mobility were constant.

Figure 3: Streamlines and decay of the current induced by a localized pump in
$d=2$ (a) below the critical point with $T=0.9T_{c}$, (b) at criticality with
$T=T_{c}$, and (c) above the critical point with $T=2T_{c}$. (d) Comparison of
the decay of the current amplitude in the direction of the pump below the
critical point (blue), at criticality (purple) and above the critical point
(red) with the expected $r^{-2}$ power law.

##  III Validity of the effective equilibrium approach

To check the validity of the effective equilibrium approach, we carried out
numerical simulations on small systems, comparing the density profile obtained
using the Kawasaki and Wolff methods (see Sec. I of the SM). The force
strength $f$ needed for the Wolff dynamics was chosen separately for each
temperature by fitting the current measurements of the Kawasaki dynamics,
shown in Fig. 3(d). As is evident from the density profiles given in Figs.
4-5, the comparison shows good agreements across the temperature range. For
larger systems with non-zero magnetization below criticality, we expect only
qualitative agreement between the two methods, as the mapping to equilibrium
breaks down.

Figure 4: Comparison of density profiles above and at the critical point with
different algorithms. The left column (a)-(d) shows the density above the
critical point $T=2T_{c}$, whereas the right column (e)-(h) shows its behavior
at criticality $T=T_{c}$. The top row ((a),(e)) shows the density profiles
measured using the Kawasaki protocol with zero average magnetization. The
second row, ((b),(f)) shows the density profiles measured using the Wolff
protocol with no offset magnetic field $h_{0}=0$. For the system above the
critical point (b), the effective magnetic field strength was set to $f=0.2$,
while in the critical point (e) it was set to $f=1.25\cdot 10^{-1}$. The third
row ((c),(g)) compares slices of the density profiles along $\hat{x}$ for
$y=0$, using the Kawasaki (blue) and Wolff (orange) protocols. (c) compares
the density profiles shown in (a) and (b), while (g) compares (e) with (f).
The bottom row, ((d),(h)) compares the angular dependence of the densities
measured on the red dashed circles appearing in the first two rows, using the
Kawasaki (blue) and Wolff (orange) protocols. (d) compares (a) with (b),
measured on a circle of radius $R=0.35L$, while (h) compares (e) with (f),
measured on a circle of radius $R=0.3L$. In both cases $L=32$. Figure 5:
Comparison of the domain walls created in the phase-separated system, below
the critical point $T=0.9T_{c}$. The domain walls predicted by the theory for
a periodic lattice, Eq. (8) and (IV), are displayed in red. The measured
domain walls appear in purple. (a) Density profile using the Kawasaki protocol
for an overall magnetization $\langle\phi\rangle=0.4$. (b) Density profile
using the Wolff protocol using the effective magnetic field of strength
$f=0.08$ and an offset magnetic field $h_{0}=1.6\cdot 10^{-3}f/a$, with $a$
the lattice constant. (c) Density profile using the Kawasaki protocol for an
overall magnetization $\langle\phi\rangle=0.2$. (b) Density profile using the
Wolff protocol using the effective magnetic field of strength $f=0.08$ and an
offset magnetic field $h_{0}=4.7\cdot 10^{-3}f/a$, with $a$ the lattice
constant. Note that despite the magnitude of the average _local_ magnetization
in each phase being very large, the equivalence between the Kawasaki and Wolff
dynamics seems to hold. The agreement is then better than expected.

##  IV Effective magnetic field for finite periodic systems

The effective magnetic field discussed in the paper was derived for an
infinite system. In two dimensions, when the pump is oriented in the $\hat{x}$
direction ${\bf f}=\hat{x}f$, the field takes the form

$h_{\rm eff}(x,y)=\frac{f}{2\pi}\frac{x}{x^{2}+y^{2}}\;.$ (2)

We remind the reader that this expression is merely the inner product of the
gradient of the Green’s function of the Laplacian and the force ${\bf f}$
exerted by the pump.

To obtain the magnetic field in a finite square periodic system of dimensions
$L\times L$, one needs to derive the Green’s function of the Laplacian in this
domain. This can be done in two steps: first, one can obtain the Green’s
function of a dipole at the origin, pointing in the $\bm{\hat{x}}$ direction
between two _infinite_ walls situated at $x=\pm L/2$. This readily gives

$h_{\parallel}(x,y;L)=\frac{f\sin\left(\frac{2\pi
x}{L}\right)}{\cosh\left(\frac{2\pi y}{L}\right)-\cos\left(\frac{2\pi
x}{L}\right)}\;,$ (3)

ensuring periodic boundary conditions along $\bm{\hat{x}}$. Second, using the
expression for the magnetic field between the infinite walls, one can obtain
the magnetic field on the square domain. To this end, we construct an infinite
array of pumps located at $(x,y)=(0,Ln)$, with $n\in\mathbb{Z}$. Using
Mathematica we obtained the following expression for the density profile:

$\displaystyle h_{\rm eff}(x,y;L)=$
$\displaystyle-i\frac{f}{L}\left[\frac{8\pi
x}{L}+\psi_{q}\left(1+i\frac{x-iy}{L}\right)\right.$
$\displaystyle-i\psi_{q}\left(1-i\frac{x+iy}{L}\right)+i\psi_{q}\left(i\frac{x+iy}{L}\right)$
$\displaystyle\left.-i\psi_{q}\left(-i\frac{x-iy}{L}\right)\right]\;$ (4)

where $\psi_{q}(z)$ is the $q$-digamma function defined by

$\psi_{q}(z)\equiv-\log(1-q)+\log(q)\sum_{n=0}^{\infty}\frac{q^{n+z}}{1-q^{n+z}}\;,$
(5)

and $q=e^{2\pi}$. We stress here that the above expression is purely real.

##  V The behavior just above the critical point

In the analysis carried out in the main text we consider a pump exactly at the
critical point where the correlation length $\xi_{0}$, in the absence of the
pump, is infinite. In this Section we discuss the behavior near the critical
point and for $T\gtrsim T_{c}$. We start by considering dimensions $d<4$.

Figure 6: The different length scales $\xi_{0}$ (orange), $\xi(h_{\rm eff})$
(red), and $r$ (blue) the distance from the pump as function of $r$ in the two
scenarios described in the text (a) $\xi_{0}>r^{*}$, (b) $\xi_{0}<r^{*}$. The
correlation length that controls the behavior of the system is marked by the
black dotted line. The regions marked in shades of gray describe the behavior
of the system – ‘NF’ for critical near-field, ‘FF’ for critical far-field, and
‘OC’ for off-critical high-temperature behavior. The red curve $\xi(h)$ is
drawn for $d=2$, having the same qualitative behavior as that of $d=3$.

Off criticality the thermal correlation length $\xi_{0}\sim\tau^{-\nu}$ is
finite. We expect that on large length-scales the $\tau>0$ behavior of the
main text is recovered. More precisely, $\phi^{\rm OC}({\rm r})\sim\chi h_{\rm
eff}({\rm r})$ where $\chi$ is the magnetic susceptibility which was given in
the main text by $\tau^{-1}$ but close to the critical point is
$\chi\sim\tau^{-\gamma}$ where $\gamma=7/4,\,1.24$ in two and three
dimensions, respectively. However, one also expects a regime where a behavior
similar to the critical behavior of the main text applies.

Recall that we showed that at $T=T_{c}$ the correlation length is finite as a
result of the magnetic field and it depends on the distance $r$ from the pump.
The interplay between the correlation length and the rate of change of $h_{\rm
eff}({\bf r})$ dictated whether the far-field or near-field behavior are
exhibited. In particular, when the correlation length at a distance $r$ from
the pump is larger than $r$ the density is controlled by the near-field
behavior.

To understand the behavior of the order parameter on different length scales,
we first have to identify the correlation length of the system as a function
of $r$. This depends on $h_{\rm eff}({\bf r})$ as well as the deviation of the
temperature from the critical value, $\tau$. As a qualitative approximation,
it is the minimum of the correlation length induced by either the magnetic
field, $\xi(h_{\rm eff}({\bf r}))$, or by a nonzero $\tau$, $\xi_{0}$.

To proceed we plot in Fig. 6 $\xi_{0},\,\xi(h_{\rm eff})$ and $r$ as functions
of $r$ using $\xi(h_{\rm eff})\sim h_{\rm
eff}^{-\frac{(\delta-1)}{\delta(2-\eta)}}$. The three functions intersect at
three points. As in the main text we denote by $r^{*}$ the intersection of $r$
with $\xi(h_{\rm eff})$. Another important intersection point, denoted by
$\overline{r}$, occurs when $\xi(h_{\rm eff})=\xi_{0}$.

Consider first the case depicted in Fig. 6(a) where $r^{*}<\overline{r}$. This
occurs when the pump is strong or $\tau$ is very small, specifically,

$\tau^{\nu(1-\frac{d}{\delta+1})}<f\;.$ (6)

In this case when $r<\overline{r}$ the correlation length which controls the
system is $\xi(h_{\rm eff})$. Then in this region the system behaves like a
critical system with a crossover from the near-field to the far field behavior
at $r^{*}$. For distances $r>\overline{r}$ the correlation length $\xi_{0}$
controls the fluctuations, resulting in a behavior similar to that of the
$\tau>0$ regime of the main text. In particular, we have $\phi^{\rm
OC}\sim\chi h_{\rm eff}({\bf r})$, as mentioned above. Indeed, it is easy to
check the self-consistency of this argument by verifying that $\phi^{\rm
OC}\sim\phi^{\rm nonlinear}$ at $r=\overline{r}$ using the scaling relation
$\gamma=\nu(2-\eta)$ [2].

Next, consider the case depicted in Fig. 6(b) where $r^{*}>\overline{r}$,
opposite to the condition in Eq. (6). In this case the correlation length on a
scale $r$ is given by $\xi(h_{\rm eff})$ for $r<\overline{r}$ and by $\xi_{0}$
for $r>\overline{r}$. Note that this implies that for $r<\xi_{0}$ the
correlation length at a distance $r$ from the pump is larger than $r$ so that
the near-field behavior is observed. For $r>\xi_{0}$ again the correlation
length can be ignored and one recovers the $\tau>0$ off-critical behavior of
the main-text. The self-consistency of the argument can by checked by
verifying that $\phi^{\rm OC}\sim\phi^{\rm linear}$ at $r=\xi_{0}$.
Interestingly, when Eq. (6) is satisfied, there are three different regimes,
but only two when it is not.

For dimensions $d>4$, the situation simplifies as there are no near and far-
field regimes. Then on scales $r$ larger than $\xi_{0}$ the density profile
matches its off-critical high-temperature form $\phi\sim{\bf f}\cdot{\bf
r}/r^{d}$ while for $r<\xi_{0}$ critical mean-field solution is valid with
$\phi\sim{\bf f}\cdot{\bf r}/r^{d-2}$.

##  VI Non-zero magnetization below $T_{c}$

Here we consider the interesting case when the average magnetization in the
system is non zero below the critical temperature. While the mapping breaks
down in this regime, a qualitative understanding can be inferred by assuming a
constant mobility. We start by considering an infinite system. Rather than
fixing the magnetization (which is not well-defined for an infinite system),
we study the shape of domains by considering the equilibrium between the two
phases. To do so we minimize the free-energy with the effective magnetic field
$h_{\rm eff}$ to get

$\displaystyle\langle\mu\rangle({\bf r})=$ $\displaystyle
K\nabla^{2}\langle\phi\rangle-\left(1+\frac{u}{\tau}\langle\phi\rangle^{2}\right)\tau\langle\phi\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{S_{d}}\frac{{\bf f}\cdot{\bf
r}}{r^{d}}-h_{0}\ .$ (7)

Here we added the field $h_{0}$ as a Lagrange multiplier which allow for a
non-zero magnetization in the system. On large length scales the gradient term
can be neglected. To identify the phase boundary, we note that the
magnetization can point in either direction when the local magnetic field,
including the $h_{0}$ is zero. Namely, when

$h_{0}=\frac{1}{S_{d}}\frac{{\bf f}\cdot{\bf r}}{r^{d}}=h_{\rm eff}\ ,$ (8)

Then the magnetization along the domain boundaries is
$\langle\phi\rangle^{2}=-\frac{\tau}{u}$. In particular, in two dimensions,
after setting ${\bf f}=\hat{x}f$, Eq. (8) reads

$\left(x-\frac{f}{4\pi h_{0}}\right)^{2}+y^{2}=\left(\frac{f}{4\pi
h_{0}}\right)^{2}\ ,$ (9)

so that a circular domain forms, passing through the pump. As in the case with
zero magnetization, the pump localizes the domain wall. For a finite size
system, when the radius of the domain reaches $\sim L/4$ this result cannot
hold anymore. However, as mentioned in Sec. IV of the SM, it is
straightforward to extend our results to account for, say, periodic boundary
conditions when deriving the expression for $h_{\rm eff}$ (IV). The results
are compared with the simulated density profiles both for Kawasaki and Wolff
simulations in Fig. 5. Note that the deviations from the theory stem from two
origins. First, the discrepancy between the Kawasaki dynamics and the Wolff
algorithm are expected because of the breaking down of the mapping to
equilibrium. The discrepancy between the theory and the Wolff algorithm is
largely due to finite-size effects and the width of the domain walls (the
agreement becomes very good for larger systems, see Fig. 7). Nonetheless, the
overall qualitative picture is correct.

Figure 7: Domain walls created in the effective equilibrium system at a
temperature below the critical point $T=0.3T_{c}$, for two values of the
offset magnetic field $h_{0}$ controlling the overall magnetization
$\langle\phi\rangle$. Density profiles, simulated using the Wolff algorithm on
a $128\times 128$ lattice, are compared with the theoretical domain walls (red
lines) predicted by Eq. (8). Note that due to the breakdown of the mapping to
effective equilibrium, these domains do not faithfully describe the original
pumped lattice gas (see Fig. 5). (a) Domain wall created for a magnetization
of $\langle\phi\rangle=5.7\cdot 10^{-2}$ due to an offset magnetic field of
$h_{0}=2.3\cdot 10^{-3}f/a$. (b) Domain wall created for a magnetization of
$\langle\phi\rangle=0.4$ due to an offset magnetic field of $h_{0}=2.3\cdot
10^{-4}f/a$.

## References

  * Newman and Barkema [1999] M. Newman and G. Barkema, Monte carlo methods in statistical physics chapter 1-4, New York, USA (1999).
  * Kardar [2007] M. Kardar,  _Statistical physics of fields_ (Cambridge University Press, 2007).
  * Wolff [1988] U. Wolff, Lattice field theory as a percolation process, Physical review letters 60, 1461 (1988).
  * Kent-Dobias and Sethna [2018] J. Kent-Dobias and J. P. Sethna, Cluster representations and the wolff algorithm in arbitrary external fields, Physical Review E 98, 063306 (2018).
  * Ferrenberg  _et al._ [2018] A. M. Ferrenberg, J. Xu, and D. P. Landau, Pushing the limits of monte carlo simulations for the three-dimensional ising model, Physical Review E 97, 043301 (2018).