* Font et al. (2011) Font, A. S., McCarthy, I. G., Crain, R. A., et al. 2011, MNRAS, 416, 2802, doi: 10.1111/j.1365-2966.2011.19227.x
  * Font et al. (2020) Font, A. S., McCarthy, I. G., Poole-Mckenzie, R., et al. 2020, MNRAS, 498, 1765, doi: 10.1093/mnras/staa2463
  * Freeman (1970) Freeman, K. C. 1970, ApJ, 160, 811, doi: 10.1086/150474
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41, doi: 10.1111/j.1365-2966.2005.09321.x
  * Geha et al. (2017) Geha, M., Wechsler, R. H., Mao, Y.-Y., et al. 2017, ApJ, 847, 4, doi: 10.3847/1538-4357/aa8626
  * Gezari et al. (2022) Gezari, S., Bentz, M., De, K., et al. 2022, arXiv e-prints, arXiv:2202.12311, doi: 10.48550/arXiv.2202.12311
  * Gilbert et al. (2019) Gilbert, K. M., Kirby, E. N., Escala, I., et al. 2019, ApJ, 883, 128, doi: 10.3847/1538-4357/ab3807
  * Gilbert et al. (2009) Gilbert, K. M., Guhathakurta, P., Kollipara, P., et al. 2009, ApJ, 705, 1275, doi: 10.1088/0004-637X/705/2/1275
  * Gilbert et al. (2012) Gilbert, K. M., Guhathakurta, P., Beaton, R. L., et al. 2012, ApJ, 760, 76, doi: 10.1088/0004-637X/760/1/76
  * Gilbert et al. (2014) Gilbert, K. M., Kalirai, J. S., Guhathakurta, P., et al. 2014, ApJ, 796, 76, doi: 10.1088/0004-637X/796/2/76
  * Gilhuly et al. (2022) Gilhuly, C., Merritt, A., Abraham, R., et al. 2022, ApJ, 932, 44, doi: 10.3847/1538-4357/ac6750
  * Governato et al. (2010) Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203, doi: 10.1038/nature08640
  * Grand et al. (2021) Grand, R. J. J., Marinacci, F., Pakmor, R., et al. 2021, MNRAS, 507, 4953, doi: 10.1093/mnras/stab2492
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, ApJ, 746, 125, doi: 10.1088/0004-637X/746/2/125
  * Harmsen et al. (2023) Harmsen, B., Bell, E. F., D’Souza, R., et al. 2023, MNRAS, doi: 10.1093/mnras/stad2480
  * Harmsen et al. (2017) Harmsen, B., Monachesi, A., Bell, E. F., et al. 2017, MNRAS, 466, 1491, doi: 10.1093/mnras/stw2992
  * Haywood et al. (2018) Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113, doi: 10.3847/1538-4357/aad235
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85, doi: 10.1038/s41586-018-0625-x
  * Helmi & White (1999) Helmi, A., & White, S. D. M. 1999, MNRAS, 307, 495, doi: 10.1046/j.1365-8711.1999.02616.x
  * Helmi et al. (1999) Helmi, A., White, S. D. M., de Zeeuw, P. T., & Zhao, H. 1999, Nature, 402, 53, doi: 10.1038/46980
  * Hood et al. (2018) Hood, C. E., Kannappan, S. J., Stark, D. V., et al. 2018, ApJ, 857, 144, doi: 10.3847/1538-4357/aab719
  * Horta et al. (2023) Horta, D., Cunningham, E. C., Sanderson, R. E., et al. 2023, ApJ, 943, 158, doi: 10.3847/1538-4357/acae87
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ibata et al. (2001a) Ibata, R., Irwin, M., Lewis, G., Ferguson, A. M. N., & Tanvir, N. 2001a, Nature, 412, 49, doi: 10.1038/35083506
  * Ibata et al. (2001b) Ibata, R., Lewis, G. F., Irwin, M., Totten, E., & Quinn, T. 2001b, ApJ, 551, 294, doi: 10.1086/320060
  * Ibata et al. (2007) Ibata, R., Martin, N. F., Irwin, M., et al. 2007, ApJ, 671, 1591, doi: 10.1086/522574
  * Ibata et al. (1994) Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194, doi: 10.1038/370194a0
  * Ibata et al. (2019) Ibata, R. A., Malhan, K., & Martin, N. F. 2019, ApJ, 872, 152, doi: 10.3847/1538-4357/ab0080
  * Ibata et al. (2014) Ibata, R. A., Lewis, G. F., McConnachie, A. W., et al. 2014, ApJ, 780, 128, doi: 10.1088/0004-637X/780/2/128
  * Ivezić et al. (2000) Ivezić, Ž., Goldston, J., Finlator, K., et al. 2000, AJ, 120, 963, doi: 10.1086/301455
  * Jang et al. (2020) Jang, I. S., de Jong, R. S., Holwerda, B. W., et al. 2020, A&A, 637, A8, doi: 10.1051/0004-6361/201936994
  * Johns et al. (2012) Johns, M., McCarthy, P., Raybould, K., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8444, Ground-based and Airborne Telescopes IV, ed. L. M. Stepp, R. Gilmozzi, & H. J. Hall, 84441H, doi: 10.1117/12.926716
  * Johnston (1998) Johnston, K. V. 1998, ApJ, 495, 297, doi: 10.1086/305273
  * Johnston et al. (2008) Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936, doi: 10.1086/592228
  * Johnston et al. (1996) Johnston, K. V., Hernquist, L., & Bolte, M. 1996, ApJ, 465, 278, doi: 10.1086/177418
  * Johnston et al. (2001) Johnston, K. V., Sackett, P. D., & Bullock, J. S. 2001, ApJ, 557, 137, doi: 10.1086/321644
  * Johnston et al. (1995) Johnston, K. V., Spergel, D. N., & Hernquist, L. 1995, ApJ, 451, 598, doi: 10.1086/176247
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864, doi: 10.1086/523619
  * Kalirai et al. (2006) Kalirai, J. S., Gilbert, K. M., Guhathakurta, P., et al. 2006, ApJ, 648, 389, doi: 10.1086/505697
  * Karachentsev et al. (2020) Karachentsev, I. D., Makarova, L. N., Brent Tully, R., et al. 2020, A&A, 643, A124, doi: 10.1051/0004-6361/202038928
  * Keller (2022) Keller, B. W. 2022, ApJ, 939, 4, doi: 10.3847/1538-4357/ac92fe
  * Keller et al. (2016) Keller, B. W., Wadsley, J., & Couchman, H. M. P. 2016, MNRAS, 463, 1431, doi: 10.1093/mnras/stw2029
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102, doi: 10.1088/0004-637X/779/2/102
  * Kravtsov et al. (2018) Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astronomy Letters, 44, 8, doi: 10.1134/S1063773717120015
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lancaster et al. (2022) Lancaster, L., Pearson, S., Williams, B. F., et al. 2022, AJ, 164, 142, doi: 10.3847/1538-3881/ac8a95
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193, doi: 10.48550/arXiv.1110.3193
  * Lee et al. (2015) Lee, D. M., Johnston, K. V., Sen, B., & Jessop, W. 2015, ApJ, 802, 48, doi: 10.1088/0004-637X/802/1/48
  * Li et al. (2022) Li, T. S., Ji, A. P., Pace, A. B., et al. 2022, ApJ, 928, 30, doi: 10.3847/1538-4357/ac46d3
  * Licquia & Newman (2015) Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96, doi: 10.1088/0004-637X/806/1/96
  * Lim et al. (2022) Lim, S. H., Raman, K. A., Buckley, M. R., & Shih, D. 2022, arXiv e-prints, arXiv:2211.11765, doi: 10.48550/arXiv.2211.11765
  * Lochhaas et al. (2021) Lochhaas, C., Tumlinson, J., O’Shea, B. W., et al. 2021, ApJ, 922, 121, doi: 10.3847/1538-4357/ac2496
  * Lochhaas et al. (2023) Lochhaas, C., Tumlinson, J., Peeples, M. S., et al. 2023, ApJ, 948, 43, doi: 10.3847/1538-4357/acbb06
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201, doi: 10.48550/arXiv.0912.0201
  * Mackereth & Bovy (2020) Mackereth, J. T., & Bovy, J. 2020, MNRAS, 492, 3631, doi: 10.1093/mnras/staa047
  * Majewski et al. (1996) Majewski, S. R., Munn, J. A., & Hawley, S. L. 1996, ApJ, 459, L73, doi: 10.1086/309944
  * Malhan & Ibata (2018) Malhan, K., & Ibata, R. A. 2018, MNRAS, 477, 4063, doi: 10.1093/mnras/sty912
  * Malin & Hadley (1997) Malin, D., & Hadley, B. 1997, PASA, 14, 52, doi: 10.1071/AS97052
  * Maller & Dekel (2002) Maller, A. H., & Dekel, A. 2002, MNRAS, 335, 487, doi: 10.1046/j.1365-8711.2002.05646.x
  * Martin et al. (2022) Martin, G., Bazkiaei, A. E., Spavone, M., et al. 2022, MNRAS, 513, 1459, doi: 10.1093/mnras/stac1003
  * Martínez-Delgado et al. (2010) Martínez-Delgado, D., Gabany, R. J., Crawford, K., et al. 2010, AJ, 140, 962, doi: 10.1088/0004-6256/140/4/962
  * Martínez-Delgado et al. (2023) Martínez-Delgado, D., Cooper, A. P., Román, J., et al. 2023, A&A, 671, A141, doi: 10.1051/0004-6361/202245011
  * McCarthy et al. (2012) McCarthy, I. G., Font, A. S., Crain, R. A., et al. 2012, MNRAS, 420, 2245, doi: 10.1111/j.1365-2966.2011.20189.x
  * McConnachie et al. (2009) McConnachie, A. W., Irwin, M. J., Ibata, R. A., et al. 2009, Nature, 461, 66, doi: 10.1038/nature08327
  * Merritt et al. (2020) Merritt, A., Pillepich, A., van Dokkum, P., et al. 2020, MNRAS, 495, 4570, doi: 10.1093/mnras/staa1164
  * Merritt et al. (2016) Merritt, A., van Dokkum, P., Abraham, R., & Zhang, J. 2016, ApJ, 830, 62, doi: 10.3847/0004-637X/830/2/62
  * Monachesi et al. (2016a) Monachesi, A., Bell, E. F., Radburn-Smith, D. J., et al. 2016a, MNRAS, 457, 1419, doi: 10.1093/mnras/stv2987
  * Monachesi et al. (2016b) Monachesi, A., Gómez, F. A., Grand, R. J. J., et al. 2016b, MNRAS, 459, L46, doi: 10.1093/mnrasl/slw052
  * Monachesi et al. (2013) Monachesi, A., Bell, E. F., Radburn-Smith, D. J., et al. 2013, ApJ, 766, 106, doi: 10.1088/0004-637X/766/2/106
  * Monachesi et al. (2019) Monachesi, A., Gómez, F. A., Grand, R. J. J., et al. 2019, MNRAS, 485, 2589, doi: 10.1093/mnras/stz538
  * Montes et al. (2023) Montes, M., Annibali, F., Bellazzini, M., et al. 2023, arXiv e-prints, arXiv:2306.09414, doi: 10.48550/arXiv.2306.09414
  * Mouhcine et al. (2005a) Mouhcine, M., Ferguson, H. C., Rich, R. M., Brown, T. M., & Smith, T. E. 2005a, ApJ, 633, 821, doi: 10.1086/468178
  * Mouhcine et al. (2005b) Mouhcine, M., Rich, R. M., Ferguson, H. C., Brown, T. M., & Smith, T. E. 2005b, ApJ, 633, 828, doi: 10.1086/468179
  * Munshi et al. (2013) Munshi, F., Governato, F., Brooks, A. M., et al. 2013, ApJ, 766, 56, doi: 10.1088/0004-637X/766/1/56
  * Mutlu-Pakdil et al. (2021) Mutlu-Pakdil, B., Sand, D. J., Crnojević, D., et al. 2021, ApJ, 918, 88, doi: 10.3847/1538-4357/ac0db8
  * Naidu et al. (2020) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48, doi: 10.3847/1538-4357/abaef4
  * Naidu et al. (2022) —. 2022, arXiv e-prints, arXiv:2204.09057, doi: 10.48550/arXiv.2204.09057
  * Newberg et al. (2002) Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245, doi: 10.1086/338983
  * Panithanpaisal et al. (2021) Panithanpaisal, N., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 920, 10, doi: 10.3847/1538-4357/ac1109
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Pearson et al. (2022) Pearson, S., Clark, S. E., Demirjian, A. J., et al. 2022, ApJ, 926, 166, doi: 10.3847/1538-4357/ac4496
  * Pearson et al. (2015) Pearson, S., Küpper, A. H. W., Johnston, K. V., & Price-Whelan, A. M. 2015, ApJ, 799, 28, doi: 10.1088/0004-637X/799/1/28
  * Pearson et al. (2019) Pearson, S., Starkenburg, T. K., Johnston, K. V., et al. 2019, ApJ, 883, 87, doi: 10.3847/1538-4357/ab3e06
  * Peeples et al. (2014) Peeples, M. S., Werk, J. K., Tumlinson, J., et al. 2014, ApJ, 786, 54, doi: 10.1088/0004-637X/786/1/54
  * Peeples et al. (2019) Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019, ApJ, 873, 129, doi: 10.3847/1538-4357/ab0654
  * Pillepich et al. (2015) Pillepich, A., Madau, P., & Mayer, L. 2015, ApJ, 799, 184, doi: 10.1088/0004-637X/799/2/184
  * Pillepich et al. (2014) Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014, MNRAS, 444, 237, doi: 10.1093/mnras/stu1408
  * Pohlen & Trujillo (2006) Pohlen, M., & Trujillo, I. 2006, A&A, 454, 759, doi: 10.1051/0004-6361:20064883
  * Pontzen & Tremmel (2018) Pontzen, A., & Tremmel, M. 2018, ApJS, 237, 23, doi: 10.3847/1538-4365/aac832
  * Purcell et al. (2010) Purcell, C. W., Bullock, J. S., & Kazantzidis, S. 2010, MNRAS, 404, 1711, doi: 10.1111/j.1365-2966.2010.16429.x
  * Purcell et al. (2007) Purcell, C. W., Bullock, J. S., & Zentner, A. R. 2007, ApJ, 666, 20, doi: 10.1086/519787
  * Purcell et al. (2008) —. 2008, MNRAS, 391, 550, doi: 10.1111/j.1365-2966.2008.13938.x
  * Querejeta et al. (2015) Querejeta, M., Meidt, S. E., Schinnerer, E., et al. 2015, ApJS, 219, 5, doi: 10.1088/0067-0049/219/1/5
  * Radburn-Smith et al. (2011) Radburn-Smith, D. J., de Jong, R. S., Seth, A. C., et al. 2011, ApJS, 195, 18, doi: 10.1088/0067-0049/195/2/18
  * Renda et al. (2005) Renda, A., Gibson, B. K., Mouhcine, M., et al. 2005, MNRAS, 363, L16, doi: 10.1111/j.1745-3933.2005.00075.x
  * Rey & Starkenburg (2022) Rey, M. P., & Starkenburg, T. K. 2022, MNRAS, 510, 4208, doi: 10.1093/mnras/stab3709
  * Robertson et al. (2005) Robertson, B., Bullock, J. S., Font, A. S., Johnston, K. V., & Hernquist, L. 2005, ApJ, 632, 872, doi: 10.1086/452619
  * Román et al. (2021) Román, J., Castilla, A., & Pascual-Granado, J. 2021, A&A, 656, A44, doi: 10.1051/0004-6361/202142161
  * Román & Trujillo (2018) Román, J., & Trujillo, I. 2018, Research Notes of the American Astronomical Society, 2, 144, doi: 10.3847/2515-5172/aad8b8
  * Román et al. (2020) Román, J., Trujillo, I., & Montes, M. 2020, A&A, 644, A42, doi: 10.1051/0004-6361/201936111
  * Sanderson et al. (2018) Sanderson, R. E., Garrison-Kimmel, S., Wetzel, A., et al. 2018, ApJ, 869, 12, doi: 10.3847/1538-4357/aaeb33
  * Sanderson et al. (2020) Sanderson, R. E., Wetzel, A., Loebman, S., et al. 2020, ApJS, 246, 6, doi: 10.3847/1538-4365/ab5b9d
  * Searle & Zinn (1978) Searle, L., & Zinn, R. 1978, ApJ, 225, 357, doi: 10.1086/156499
  * Sesar et al. (2011) Sesar, B., Jurić, M., & Ivezić, Ž. 2011, ApJ, 731, 4, doi: 10.1088/0004-637X/731/1/4
  * Shang et al. (1998) Shang, Z., Zheng, Z., Brinks, E., et al. 1998, ApJ, 504, L23, doi: 10.1086/311563
  * Shankar et al. (2006) Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006, ApJ, 643, 14, doi: 10.1086/502794
  * Sharma et al. (2020) Sharma, R. S., Brooks, A. M., Somerville, R. S., et al. 2020, ApJ, 897, 103, doi: 10.3847/1538-4357/ab960e
  * Shih et al. (2022) Shih, D., Buckley, M. R., Necib, L., & Tamanas, J. 2022, MNRAS, 509, 5992, doi: 10.1093/mnras/stab3372
  * Shipp et al. (2018) Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862, 114, doi: 10.3847/1538-4357/aacdab
  * Shipp et al. (2023) Shipp, N., Panithanpaisal, N., Necib, L., et al. 2023, ApJ, 949, 44, doi: 10.3847/1538-4357/acc582
  * Sick et al. (2015) Sick, J., Courteau, S., Cuillandre, J.-C., et al. 2015, in Galaxy Masses as Constraints of Formation Models, ed. M. Cappellari & S. Courteau, Vol. 311, 82–85, doi: 10.1017/S1743921315003440
  * Simons et al. (2020) Simons, R. C., Peeples, M. S., Tumlinson, J., et al. 2020, ApJ, 905, 167, doi: 10.3847/1538-4357/abc5b8
  * Simpson et al. (2013) Simpson, C. M., Bryan, G. L., Johnston, K. V., et al. 2013, MNRAS, 432, 1989, doi: 10.1093/mnras/stt474
  * Smercina et al. (2022) Smercina, A., Bell, E. F., Samuel, J., & D’Souza, R. 2022, ApJ, 930, 69, doi: 10.3847/1538-4357/ac5d56
  * Smercina et al. (2020) Smercina, A., Bell, E. F., Price, P. A., et al. 2020, ApJ, 905, 60, doi: 10.3847/1538-4357/abc485
  * Smith et al. (2011) Smith, B. D., Hallman, E. J., Shull, J. M., & O’Shea, B. W. 2011, ApJ, 731, 6, doi: 10.1088/0004-637X/731/1/6
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217, doi: 10.1093/mnras/stw3291
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757, doi: 10.48550/arXiv.1503.03757
  * Stinson et al. (2010) Stinson, G. S., Bailin, J., Couchman, H., et al. 2010, MNRAS, 408, 812, doi: 10.1111/j.1365-2966.2010.17187.x
  * Streich et al. (2014) Streich, D., de Jong, R. S., Bailin, J., et al. 2014, A&A, 563, A5, doi: 10.1051/0004-6361/201220956
  * Thomas et al. (2003) Thomas, D., Maraston, C., & Bender, R. 2003, MNRAS, 339, 897, doi: 10.1046/j.1365-8711.2003.06248.x
  * Tissera et al. (2013) Tissera, P. B., Scannapieco, C., Beers, T. C., & Carollo, D. 2013, MNRAS, 432, 3391, doi: 10.1093/mnras/stt691
  * Tollerud et al. (2011) Tollerud, E. J., Boylan-Kolchin, M., Barton, E. J., Bullock, J. S., & Trinh, C. Q. 2011, ApJ, 738, 102, doi: 10.1088/0004-637X/738/1/102
  * Trujillo & Fliri (2016) Trujillo, I., & Fliri, J. 2016, ApJ, 823, 123, doi: 10.3847/0004-637X/823/2/123
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, The Astrophysical Journal Supplement Series, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * Unavane et al. (1996) Unavane, M., Wyse, R. F. G., & Gilmore, G. 1996, MNRAS, 278, 727, doi: 10.1093/mnras/278.3.727
  * Venn et al. (2004) Venn, K. A., Irwin, M., Shetrone, M. D., et al. 2004, AJ, 128, 1177, doi: 10.1086/422734
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Weinberg et al. (2023) Weinberg, D. H., Griffith, E. J., Johnson, J. W., & Thompson, T. A. 2023, arXiv e-prints, arXiv:2309.05719.  https://arxiv.org/abs/2309.05719
  * Weisz et al. (2011) Weisz, D. R., Dalcanton, J. J., Williams, B. F., et al. 2011, ApJ, 739, 5, doi: 10.1088/0004-637X/739/1/5
  * WFIRST Astrometry Working Group et al. (2019) WFIRST Astrometry Working Group, Sanderson, R. E., Bellini, A., et al. 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 044005, doi: 10.1117/1.JATIS.5.4.044005
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Worthey (1994) Worthey, G. 1994, ApJS, 95, 107, doi: 10.1086/192096
  * Xue et al. (2015) Xue, X.-X., Rix, H.-W., Ma, Z., et al. 2015, ApJ, 809, 144, doi: 10.1088/0004-637X/809/2/144
  * Yanny et al. (2000) Yanny, B., Newberg, H. J., Kent, S., et al. 2000, ApJ, 540, 825, doi: 10.1086/309386
  * Yoachim (2022) Yoachim, P. 2022, 016: Surface brightness limit derivations, LSST.  https://smtn-016.lsst.io/#change-record
  * Yu et al. (2020) Yu, S., Bullock, J. S., Wetzel, A., et al. 2020, MNRAS, 494, 1539, doi: 10.1093/mnras/staa522
  * Zheng et al. (2020) Zheng, Y., Peeples, M. S., O’Shea, B. W., et al. 2020, ApJ, 896, 143, doi: 10.3847/1538-4357/ab960a
  * Zolotov et al. (2009) Zolotov, A., Willman, B., Brooks, A. M., et al. 2009, ApJ, 702, 1058, doi: 10.1088/0004-637X/702/2/1058


# Simulated surface diffusion in nanoporous gold and its dependence on surface
curvature

Conner Marie Winkeljohn<EMAIL_ADDRESS>Sadi Md Shahriar  Erkin
Seker  Jeremy K. Mason<EMAIL_ADDRESS>Department of Materials Science
and Engineering, University of California, Davis, Davis, CA 95616, USA
Department of Electrical and Computer Engineering, University of California,
Davis, Davis, CA 95616, USA

###### Abstract

The morphological evolution of nanoporous gold is generally believed to be
governed by surface diffusion. This work specifically explores the dependence
of mass transport by surface diffusion on the curvature of a gold surface. The
surface diffusivity is estimated by molecular dynamics simulations for a
variety of surfaces of constant mean curvature, eliminating any chemical
potential gradients and allowing the possible dependence of the surface
diffusivity on mean curvature to be isolated. The apparent surface diffusivity
is found to have an activation energy of $\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.74$ eV with a weak dependence on
curvature, but is consistent with the values reported in the literature. The
apparent concentration of mobile surface atoms is found to be highly variable,
having an Arrhenius dependence on temperature with an activation energy that
also has a weak curvature dependence. These activation energies depend on
curvature in such a way that the rate of mass transport by surface diffusion
is nearly independent of curvature, but with a higher activation energy of
$\raise 0.73193pt\hbox{$\scriptstyle\sim$}1.01$ eV. The curvature dependencies
of the apparent surface diffusivity and concentration of mobile surface atoms
is believed to be related to the expected lifetime of a mobile surface atom,
and has the practical consequence that a simulation study that does not
account for this finite lifetime could underestimate the activation energy for
mass transport via surface diffusion by $\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.27$ eV.

###### keywords:

nanoporous gold , surface diffusion , molecular dynamics

##  1 Introduction

Nanoporous gold (np-Au) has a bicontinuous structure with interconnected pores
and ligaments that have mean widths from nanometers to microns [1]. The
interface between the ligaments and pores create a large surface area, and
given the characteristic length scale, large surface curvatures as well. These
unique features give rise to many of the desirable properties of np-Au,
resulting in its use in a variety of applications including catalysis [2, 3,
4], optical sensors [5, 6, 7], biomedical devices [8, 9, 10], and battery
electrodes [11, 4, 12]. One reason for this diversity of uses is that many of
its notable properties (e.g., optical properties and mechanical strength) are
strongly dependent on a microstructure that can be tailored to each
application. Specifically, the means and standard deviations of the
distributions of pore and ligament widths, and the anisotropy of all four of
these quantities, can be used to modify the properties of np-Au. The
implication is that precisely controlling the morphology and topology of np-Au
is of considerable practical interest.

np-Au is produced by chemically or electrochemically dealloying a silver-gold
alloy. There has already been extensive research on controlling the np-Au
morphology by means of the fabrication parameters, including the choice of
etching solution and the initial precursor alloy constituents [13, 14, 15, 16,
2]. While the pore-to-ligament size ratio is an essential microstructural
feature, it is highly sensitive to these parameters and therefore difficult to
predict [17]. This has encouraged the community to develop ways to adjust the
pore and ligament sizes after fabrication by means of thermally-driven
coarsening, which is the most commonly used technique for this purpose.

The vast majority of studies quantify np-Au microstructure by means of SEM
images [18, 19, 1, 20] from which the width and length distributions of the
ligaments and pores can be calculated. One method of directly calculating
ligament widths and lengths involves reducing the structure pixel by pixel to
a skeletal system, counting the number of deleted pixels between the surface
and the skeleton to determine ligament radius, and then measuring the length
of the skeleton branches [17, 21]. An alternative is to take the length of the
thinnest part of each ligament as the ligament width [4]. The literature on
this subject is complicated not only by this diversity of definitions, but by
the infrequent reporting of standard deviations of ligament width
distributions [1].

A small number of studies use 3D imaging techniques to quantify other aspects
of the np-Au microstructure including the pore and ligament topology and the
distribution of surface curvatures. The surface curvature at any point can be
described by the two principal curvatures $\kappa_{1}$ and $\kappa_{2}$ which
are defined as the inverse radii of curvature along the two principal
directions, with $\kappa_{1}$ as the principal curvature of greatest
magnitude. Chen et al. [22] used X-ray tomography to image structures in 3D
and generate a surface mesh, allowing them to find the surface area per volume
$S_{A}$ and the principal curvatures at each surface point. The probability
distribution of finding a pair of principal curvatures at a surface point is
called the interface shape distribution (ISD) and is perhaps the standard way
to analyze such curvature information [22, 23, 18]. The ISD appears to be
correlated with other geometric quantities, with $S_{A}$ generally increasing
and the volume fraction of pores generally decreasing with the surface mean
curvature [22, 17]. Part of the reason such a wide variety of characterization
measurements are used is that the key characteristics which influence
morphology evolution have yet to be precisely identified and related to easily
measurable features.

The kinetics of thermally-driven coarsening are usually measured by analyzing
SEM images before and after thermal annealing. The consensus in the literature
is that this coarsening occurs by surface diffusion, with prior investigations
having ruled out bulk diffusion as well as evaporation/condensation mechanisms
[24, 16, 22, 25, 21]. The mass transport is driven by a chemical potential
gradient that results from variations in the surface curvature as given by the
Gibbs–Thompson equation; for a surface of mean curvature
$H=\kappa_{1}+\kappa_{2}$, the chemical potential $\mu(H)$ of an atom on the
curved surface can be related to that of an atom on a flat surface $\mu(0)$ by

$\mu(H)=\mu(0)+\gamma\Omega H$ (1)

where $\gamma$ is the surface energy per area and $\Omega$ is the atomic
volume. Kolluri and Demkowicz [26] instead suggest that coarsening could occur
by means of network restructuring where localized plasticity allows ligaments
to collapse onto their neighbors, making the rate of coarsening dependent on
the relative density. However, they considered a structure with grain sizes
comparable to the ligament width, whereas others have observed much larger
grain sizes in np-Au [15].

There is a frequently-investigated hypothesis that thermal coarsening of np-Au
leads to statistically self-similar microstructures [27, 26, 23, 28, 21].
Given that surface diffusion is the primary mass transport mechanism, a
scaling argument suggests that every characteristic length scale $\lambda$ of
a statistically self-similar microstructure increases with time $t$ as:

$\lambda\propto\big{(}tD_{S}\big{)}^{n}$ (2)

where $n$ is the coarsening exponent and $D_{S}$ is the surface diffusivity
[29, 21]. Measuring $D_{S}=D_{0}\exp[-E_{a}/(k_{B}T)]$ as a function of
temperature would then allow the activation energy $E_{a}$ for coarsening to
be found. If this relationship holds, the characteristic length scale of np-Au
could be easily and predictably controlled by changing the temperature or
duration of a thermal anneal.

Efforts in the literature along these lines have yielded conflicting results
though, with no consensus that thermal coarsening is self similar and a wide
variety of reported activation energies. Early kinetic Monte Carlo (kMC)
simulations found that self-similar evolution was unlikely ($E_{a}=0.9$ eV)
[30], but recent kMC simulations concluded that coarsening was self similar
($E_{a}=0.65$ eV) based on an analysis of ligament spacing [25]. An early
experimental study found that the initial coarsening was not self-similar
($E_{a}=0.643$ eV), [23], but Jeon et al. concluded that coarsening was self
similar ($E_{a}=0.353$ eV) based on 3D reconstructions of np-Au
microstructures [31]. Whereas the above studies used coarsening exponents of
$n=4$, McCue et al. reported self-similar coarsening with $n=8$ based on an
analysis of published 2D images [1], and Gao et al. instead reported
coarsening that was not self similar with $n=3.5$ [21]. Finally, a study by
Ziehmer et al. using 3D reconstructions to analyze both principal curvatures
was inconclusive as to whether the coarsening of np-Au could be described as
self-similar, stating that different regimes could characterize early and late
stage coarsening [18].

Our conclusion is that thermal coarsening is still likely to be governed by
surface diffusion, but that the factors influencing surface diffusion are not
sufficiently well understood for models to accurately predict experimental
coarsening rates. For example, if the surface diffusivity depends on the
surface curvature, then $D_{S}$ in Eq. 2 would be a function of $\lambda$.
Equation 2 whould then not be expected to hold for any fixed value of $D_{S}$
during an interval where $\lambda$ changes appreciably, perhaps leading to the
conclusion that thermal coarsening was not self-similar. Such a possibility
would be most likely to occur when $\lambda$ is small since it is for those
values that the surface curvature most rapidly changes with respect to
$\lambda$. Many of the conflicting results above could then be explained as
consequences of differing initial ligament width distributions.

Our intention is to investigate whether the surface diffusivity could in fact
depend on the mean curvature of the surface, at least for the ligaments and
pores in np-Au. While there have been several kMC studies considering the
coarsening of np-Au [30, 25], these techniques require that the relevant
atomic migration rates be specified beforehand. We instead use molecular
dynamics (MD) to study surface diffusion while making minimal prior
assumptions about the relevant atomic population and migration mechanisms.
Section 2 describes both the setup and conduct of the MD simulations and our
approach to extracting surface diffusion coefficients. The results of this
analysis are discussed in Sec. 3, where the surface diffusivity and
concentration of mobile surface atoms appear to have a slight dependence on
the surface mean curvature. This is believed to be related to the mobile
surface atoms having a finite curvature-dependent lifetime [32] that is not
accounted for by the analysis, with the true concentration of mobile surface
atoms being constant and independent of curvature. A finite mobile surface
atom lifetime would increase the apparent surface atom population (the
population of atoms that diffuse on the surface at any point during the
simulation) and decrease the apparent surface diffusivity (a mobile surface
atom migrates for only a fraction of the simulation time) compared to the true
values of these quantities. The lifetime is expected to decrease with an
increasing rate of exchange between the mobile surface atom population and a
second population of comparatively immobile atoms, e.g., ones in the bulk. A
curvature-independent rate of mass transport by surface diffusion then implies
that the true surface diffusivity is independent of curvature as well. This
raises the general concern that an analysis that does not account for the
finite curvature-dependent lifetimes of the mobile surface atoms could find an
activation energy for mass transport that differs by $\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.27$ eV from the values reported here. The
consequences of this for predicting experimental np-Au coarsening rates are
discussed with the conclusions in Sec. 4.

##  2 Methods

The purpose of this work is to characterize diffusion on the surface of the
ligaments in np-Au, with the intention of eventually using this knowledge to
develop predictive models for thermal coarsening and general morphology
evolution. There are several reasons why measuring surface diffusivities by
means of MD simulations is difficult though. First, the requirement for
numerical stability limits the time step to about a femtosecond and the
overall length of an MD simulation to a few hundreds of nanoseconds, whereas
diffusion experimentally occurs on the order of seconds. While surface
diffusion can be accelerated by increasing the temperature, this can be done
only up to the melting point of the material without changing the mass
transport mechanism. Second, diffusion is by definition a continuum-level
phenomenon, involving atomic concentrations and chemical potentials that are
not easily defined at the atomic scale. Measuring surface diffusion rates in
MD therefore requires a sufficient number of atomic migration events to
compare the resulting probability distribution of displacements to that
predicted by the phenomenological diffusion equations. This is difficult to
achieve in simulations that are often limited by the available computational
resources to fewer than a million atoms. Third, mass transport can be driven
by gradients in any one of a variety of intensive thermodynamic quantities,
not just by a chemical potential gradient. If such effects are not accounted
for and controlled in the simulation design, the analysis of the MD results
could be subject to systematic errors that would confound subsequent modeling
efforts.

###  2.1 Molecular dynamics simulations

Figure 1: The simulation cell geometries used to extract surface diffusivity
on concave ($\kappa_{1}<0$) and convex ($\kappa_{1}>0$) surfaces. $\kappa_{2}$
is along $\hat{\boldsymbol{z}}$ and $\kappa_{1}$ along
$\hat{\boldsymbol{\theta}}$. The flat plate used as a zero mean curvature
surface is not shown.

With a view to the considerations above, this study uses three simulation cell
geometries with positive, zero, and negative surface curvatures for our MD
simulations of surface diffusion. Since the simulated atomic displacements
will eventually be compared with the predictions of the diffusion equation, it
is desirable for the analytical predictions to be as simple as possible. In
the absence of drift, diffusion theory predicts that the distribution of
atomic displacements should be identical to those for a random walk. Hence, it
is desirable that the diffusion potential should be constant everywhere, and
this effectively requires that the surface be one of constant mean curvature.
Given the characteristic ligaments and pores of the physical system, this
suggests that the natural simulation cell geometries are cylinders of positive
or negative curvature with $\kappa_{2}=0$ and $\kappa_{1}\neq 0$ as in Fig. 1.
The available computational resources limited the maximum size of the
simulation cell and therefore the cylinder radius, resulting in comparatively
high mean curvature magnitudes except for one simulation cell containing a
flat slab of Au for which the mean curvature was zero. The five simulation
cell geometries that were considered included punch holes of radius $30\
\text{\AA}$ and $34\ \text{\AA}$, a flat slab $35\ \text{\AA}$ thick, and
cylinders of radius $34\ \text{\AA}$ and $30\ \text{\AA}$, roughly
corresponding to the lower limit of feature sizes experimentally observed in
np-Au. These had corresponding nonzero principal curvatures of $-0.033$,
$-0.029$, $0.0$, $0.029$, and $0.033\ \text{\AA}^{-1}$.

The $z$-axis for the simulation box was always parallel to the axis of the
punch holes and cylinders. This had the advantage that the directions of
principal curvature in all simulations were aligned with the
$\hat{\boldsymbol{z}}$ and $\hat{\boldsymbol{\theta}}$ directions, and that
the distribution of atomic displacements could be decomposed as a product of
independent distributions along these two directions. The diffusivity is
calculated in both the $\hat{\boldsymbol{z}}$ and $\hat{\boldsymbol{\theta}}$
directions from these respective displacement distribution data sets for each
curvature.

Figure 2: The four different crystal orientations used in our simulations. The
respective crystallographic directions were oriented along the $x$-, $y$-, and
$z$-axes of the simulation cell (the cylinder axis was always along the
$z$-axis). For example, the [$211$] [$\bar{1}42$] [$\bar{2}13$] orientation
aligns the crystallographic direction [$211$] with the $x$-axis of the
simulation cell, [$\bar{1}42$] with the $y$-axis, and [$\bar{2}13$] with the
$z$-axis.

The diffusivity is modeled as being independent of crystallographic
orientation given that such a quantity is more likely to be relevant to the
coarsening of np-Au when averaged over crystal orientations. Of course, the
crystallographic orientation and the coordination of surface sites does have a
significant impact on the surface diffusivity that needs to be accounted for
[33, 34]. Simulating just a single crystallographic orientation would fix the
surface facets and potentially bias the coordination distribution’s dependence
on curvature and the apparent surface diffusivity. That is, simply varying the
radius of a pillar with a fixed crystallographic orientation would result in
behavior that depends on the specific choice of orientation and the facet
pattern rather than on $\kappa_{1}$ alone. Figure 2 shows that cylinders with
different crystallographic orientations can have dramatically different
patterns of surface facets and ridges, with the potential energy of the atom
and therefore the activation energy for atomic jumping depending strongly on
the type of site. Moreover, the size of atomic jumps along the
$\hat{\boldsymbol{z}}$ and $\hat{\boldsymbol{\theta}}$ directions depends
explicitly on the corresponding crystallographic directions, causing the
distribution of atomic displacements to be discrete rather than continuous as
would be predicted by continuum diffusion theory. Establishing the effective
surface diffusivity therefore requires averaging over several crystal
orientations.

The periodic boundary conditions require that the simulation cell edge lengths
be integer multiples of the crystal lattice periodicity, effectively
constraining the crystallographic directions along the coordinate axes to
those with single-digit Miller indices. The four crystallographic orientations
shown in Fig. 2 were used for the cylinders, punch holes, and slabs, with the
listed crystallographic directions aligned with the $x$-, $y$-, and $z$-axes
of the simulation cell. These orientations were selected partly to give a
variety of angles between the $\langle 111\rangle$ directions and the
$z$-axis, since diffusive jumps are expected to be easiest on $\\{111\\}$
planes. The displacement distributions of all four crystallographic
orientations were combined into a single data set for a given $\kappa_{1}$ to
find the orientation-averaged diffusivity in the subsequent analysis. This had
the added advantage of increasing the size of the sampled data to help
mitigate random errors. A comparison between the displacements from a single
simulation and those compiled from four simulations, each with a different
crystal orientation, is show in Fig. 3.

All simulations were run using a modified embedded atom method (MEAM)
potential in LAMMPS [35]. The MEAM potential was chosen over more commonly-
used EAM potentials for several reasons. As our study is concerned with the
behavior of surface atoms, an interatomic potential with inaccurate surface
energies could lead to surface reconstruction and a dramatic change to the
surface diffusivity. While standard EAM potentials [36] show good agreement to
bulk experimental properties, they are notorious for their poor agreement with
experimental surface energies (though recently there has been improvement
[37]). Whereas EAM potentials assume that the electron density $\rho$ is
spherically symmetric, the MEAM potential is more flexible and allows $\rho$
to depend on angular coordinates [38]. Since our simulations are not only
concerned with but dominated by the surface properties of Au rather than the
bulk properties, this angular dependence could be crucial for accurately
modeling the intended behavior. Numerical results supporting the selection of
the MEAM potential are discussed further in Sec. 3.

All simulations used periodic boundary conditions in all directions and were
initialized to be roughly $9$-$10$ nm in the $z$ direction. The punch hole
geometry simulations had $x$ and $y$ lengths of roughly $16$-$20$ nm and
contained $106,000$ to $160,000$ atoms, with the larger simulation size
required to prevent the periodic images of the punch holes from interacting
elastically. By comparison the cylinder geometry simulations contained only
$15,000$ to $20,000$ atoms and the slab geometry simulations contained
$19,000$ to $22,000$ atoms, and both of these had $x$ and $y$ lengths of
roughly $9$ nm. After being initialized, all simulations were relaxed by a
potential energy minimization in LAMMPS. This was followed by equilibration
for $1$ ns in the isothermal-isobaric ensemble at zero pressure and at the
same temperature subsequently used in the production runs to allow for surface
reconstruction. All simulations were finally run for $10$ ns with NPT
integration using the Nosé–Hoover thermostat and barostat, with the atomic
displacement data over this interval used for the analysis reported in Sec. 3.

###  2.2 Fitting diffusion coefficients

Surfaces are well known to be short circuit diffusion pathways, and the
diffusivity of atoms at the surface is expected to be much higher than the
diffusivity of atoms in the bulk. Estimating the diffusivity of the surface
atoms therefore requires a way to distinguish between these two populations,
ideally while introducing as little bias as possible. For example, it is not
sufficient to simply identify all of the atoms with low coordination numbers
at the beginning of the simulation as surface atoms since some of these could
eventually become bulk atoms, after which point their behavior should not be
included in the statistics. Nor is it sufficient to include atoms in the
statistics only during periods for which their coordination numbers are below
an arbitrarily specified cutoff, since the experimental surface diffusivity
describes the average behavior of atoms with a variety of environments.

Rather than first identifying the population of surface atoms and then
analyzing the behavior of this population, the approach developed here is to
do the opposite. The motion of an individual atom in the absence of a
diffusion potential gradient is reasonably modeled as that of a random walker
whose one-dimensional displacement $x$ after a time $t$ is distributed as

$P(x,t)=\frac{1}{2\sqrt{\pi Dt}}\exp\left(-\frac{x^{2}}{4Dt}\right)$ (3)

where $D$ is the temperature-dependent diffusivity. Assume that the sets of
atoms belonging to the bulk and surface populations are unknown but fixed
(this will be reconsidered below), let $D_{S}$ and $D_{B}$ respectively be the
diffusivities of the surface and bulk populations, and let $n_{S}$ be the
unknown fraction of atoms that are included in the surface population. Then
the probability distribution for the one-dimensional displacement $x$ of an
atom selected uniformly at random from all atoms in the simulation cell is

$\displaystyle P(x,t)$ $\displaystyle=\frac{n_{S}}{2\sqrt{\pi
D_{S}t}}\exp\left(-\frac{x^{2}}{4D_{S}t}\right)$
$\displaystyle\quad+\frac{1-n_{S}}{2\sqrt{\pi
D_{B}t}}\exp\left(-\frac{x^{2}}{4D_{B}t}\right).$ (4)

Observe that the quantities $D_{S}$, $D_{B}$, and $n_{S}$ can be determined by
fitting this equation to the distribution of all atomic displacements without
any of the complications described above. Indeed, since diffusion is a
continuum-level phenomenon, one expects to be able to measure surface
diffusivities without needing to specify which atoms during which periods of
time are considered as belonging to the surface population.

Figure 3:  Distributions of atomic displacements in the $z$ direction after
$10\ \mathrm{ns}$ at $900\ \mathrm{K}$ for cylinders of radius $34\
\text{\AA}$; the vertical red lines indicate the interatomic distance.
Displacements on a cylinder with [$001$] along the $z$ direction (top)
generally occur at integer multiples of a lattice translation vector.
Combining the displacement data from simulations using the four orientations
shown in Fig. 2 reduces the strength of these peaks. Figure 4: The $z$
components (top) and magnitudes (bottom) of atomic displacements on cylinders
of radius $30\ \text{\AA}$ during a $10$ ns simulation at $900$ K. The
displacement magnitude and the number of displaced atoms visibly depends on
the crystallographic orientation (left and right).

Of course, there are limitations to this approach. Equation 2.2 should only be
used to analyze displacements along an eigendirection of the surface
diffusivity tensor, since it is only under these circumstances that the
overall displacement of a random walker can be decomposed into a product of
distributions in the form of Eq. 3. The surface diffusivity is expected to
have eigendirections aligned with the directions of the principle curvatures,
though this is a reasonable assumption only if the surface diffusivity is
defined as an average over crystal orientations. Observe that in the top of
Fig. 3, the distribution of $z$-displacements has a series of peaks
corresponding to discrete jumps between atomic sites along the [$001$]
direction; the orientation of the crystal lattice generally has a strong
effect on the distribution of atomic displacements in any particular MD
simulation, emphasizing the need to average over multiple crystal
orientations. This is more directly visible in Fig. 4, where the atomic
displacements visibly depend on the crystallographic orientation (all other
simulation properties are the same). Averaging over crystal orientations as in
the bottom of Fig. 3 not only mitigates this effect, but makes clear that the
probability distribution can reasonably be modeled as a sum of Gaussians, with
a narrow Gaussian (the peak is above the vertical limit of the image)
corresponding to the bulk population and a wider Gaussian corresponding to the
surface population. For emphasis, fitting this probability distribution to Eq.
2.2 allows $D_{S}$, $D_{B}$, and $n_{S}$ to be found without needing to
provide an explicit atomic-level prescription for what constitutes a surface
atom.

Notice that a finite mobile surface atom lifetime would cause the apparent
values of $n_{S}$ and $D_{S}$ to depend on the simulation duration, with
$n_{S}$ increasing (as atoms join and leave the mobile surface atom
population) and $D_{S}$ decreasing (as atoms spend a decreasing fraction of
the simulation in the mobile state) with increasing simulation duration.
Increasing the simulation size and decreasing the duration to measure the
instantaneous values of $n_{S}$ and $D_{S}$ is not feasible since the
diffusion length needs to be at least several interatomic spacings for the
distribution of atomic displacements to be reasonably approximated as a
Gaussian distribution. That said, the time-averaged flux of surface atoms over
the simulation duration is equal to the average of the instantaneous fluxes,
and since the instantaneous fluxes should be constant by time translation
symmetry, the average flux should not depend on the simulation duration. The
average flux is also proportional to the product of the average concentration
and average velocity of the diffusing species, meaning that the product
$n_{S}D_{S}$ is expected to be independent of simulation duration even when
$n_{S}$ and $D_{S}$ individually are not. This idea is explicitly tested in
Sec. 3.2 below, and supports the existence of a finite mobile surface atom
lifetime [32].

Equation 2.2 was fit to the simulation data by maximizing the log likelihood
function

$\ell(D_{S},D_{B},n_{S})=\sum_{i}\log P(x_{i},t|D_{S},D_{B},n_{S})$ (5)

with respect to the parameters $D_{S}$, $D_{B}$, and $n_{S}$, where $x_{i}$ is
the displacement of the $i$th atom in either the $\hat{\boldsymbol{z}}$ or
$\hat{\boldsymbol{\theta}}$ direction after $10\ \mathrm{ns}$. This equation
can be explicitly differentiated with respect to these quantities, allowing
the maximization to be performed using a modified Brent’s method [39]. The
uncertainty in the parameter estimates was found using bootstrapping where an
auxiliary data set is generated by sampling from the original data set with
replacement and the parameters are estimated using the auxiliary data set.
This procedure is repeated a sufficient number of times to be able to
characterize the mean and variance of the resulting distributions of estimated
parameter values. This was implemented using $80$ auxiliary data sets, each of
the same size as the original data set, and the resulting distributions of
$D_{S}$, $D_{B}$, and $n_{S}$ are reported in Sec. 3 below.

###  2.3 Activation Energy and Temperature

The values of $D_{S}$ and $n_{S}$ were found by the procedure in Sec. 2.2 for
a variety of temperatures. $D_{S}$ is modeled as having a standard Arrhenius
dependence on temperature:

$D_{S}=D_{S}^{0}\exp\left(-\frac{\Delta H_{S}}{k_{B}T}\right)$ (6)

where $D_{S}^{0}$ is the diffusion prefactor, $\Delta H_{S}$ is the activation
enthalpy, $k_{B}$ is Boltzmann’s constant, and $T$ is the absolute
temperature. Notice that $D_{S}^{0}$ depends on a distribution of surface atom
configurations and migration events, and since these change with the surface
mean curvature, $D_{S}^{0}$ could depend on the surface mean curvature as
well. The fraction of surface atoms $n_{S}$ does not have a standard form but
is modeled by:

$n_{S}=\alpha_{0}\frac{2r_{a}A}{V}\exp\left(-\frac{E_{F}}{k_{B}T}\right)$ (7)

where $E_{F}$ is the formation energy of a mobile surface atom, $A$ and $V$
are the total surface area and volume of the gold structure (e.g., cylinder),
respectively, and $r_{a}$ is the radius of a gold atom. The prefactor can be
motivated by observing that $2r_{a}A/V$ is roughly the fraction of atoms in
the simulation cell that belong to a single surface monolayer. This implies
that $\alpha_{0}$ is a dimensionless factor proportional to the number of
surface monolayers from which an atom could become a mobile surface atom.

The activation energies and prefactors in Eqs. 6 and 7 can be found from the
slopes and intercepts of semi-log plots of $D_{S}$ and $n_{S}$ as functions of
temperature. Simulations were performed for each combination of curvature and
crystal orientation at four different temperatures, namely, $900$, $950$,
$1000$, and $1050\ \mathrm{K}$. This temperature range was selected for being
the highest possible without initiating surface melting and changing the mass
transport mechanism, the intention being to increase the number of atomic
jumps observed within a fixed simulation duration of $10\ \mathrm{ns}$. For
comparison, coarsening of experimental np-Au samples generally occurs at
temperatures of $500$-$700\ \mathrm{K}$ over times on the order of minutes to
hours [19].

##  3 Results and Discussion

###  3.1 Atomic Potentials

As stated earlier, the surface energy of np-Au could substantially affect the
mass transport mechanism. This requires careful selection of the interatomic
potential since one where the surface energy is unrealistically large could
drive the generation of dislocations and mass transport by plastic
deformation. Alternatively, an interatomic potential where the binding energy
of surface atoms is unrealistically small could increase the surface atom
mobility as a consequence of premature surface melting.

Figure 5: Atomic displacements during $1$ ns as simulated using the EAM
potential (left) or the MEAM potential (right). Both simulations started with
identical atomic coordinates for a cylinder of radius $30\ \text{\AA}$ and a
[$1\bar{7}4$] [$\bar{1}12$] [$311$] orientation and were relaxed for $10$ ns
at $1000$ K. Figure 6:  The Au cylinder on the left shows the atomic
coordination for each atom after a $10$ ns simulation at $1000$ K using the
EAM potential (left) or the MEAM potential (right).

Two otherwise identical simulations were performed to compare the effects of
using an EAM potential [36] and a MEAM potential [38] for gold. Both
simulations contained a pillar of $13,853$ gold atoms, with the pillar axis
along the $z$ axis, orientated such that the crystallographic directions
[$1\bar{7}4$], [$\bar{1}12$], and [$311$] were respectively aligned with the
$x$, $y$, and $z$ axes of the simulation box. Both simulations were
thermostated at $1000$ K and run for a total of $11$ ns, with the atomic
displacements over the final $1$ ns of the simulation plotted in Fig. 5. It is
noteworthy that the atoms in the simulation using the EAM potential can be
seen by the common neighbor analysis in Fig. 6 to be significantly more
disordered further from the surface compared to those in the simulation using
the MEAM potential. Along with the observation that atoms displace in the EAM
simulation significantly deeper from the surface, this suggests that there
could be surface melting at temperatures below the bulk melting point of gold
in the EAM simulations.

Figure 7: Histogram showing the distribution of atomic displacements in the
$\hat{\boldsymbol{\theta}}$ direction after $t=10$ ns, with one series of
varying curvature (left) and one of varying temperature (right). The
distribution of surface atom displacements measurably broadens with increasing
curvature and temperature (up the figure). The figures have been cropped in
the vertical direction since the number of negligible displacements of bulk
atoms is two orders of magnitude higher.

The model we are using to extract diffusion coefficients from atomic
displacements requires combining data from enough simulations with different
crystallographic orientations to give a displacement distribution that
approximates a smooth Gaussian. Figure 7 shows several such distributions,
each including data from all four orientations in Fig. 2 simulated at the
indicated temperature and curvature; these four orientations appear to be
sufficient to give a distribution that resembles a Gaussian, allowing Eq. 2.2
to be reliably fit to the data. The histograms are cropped to better show the
part of the distribution relevant to surface diffusion as there are orders of
magnitude more counts for the negligible displacements of the bulk atoms.
Observe that the histograms measurably broaden with increasing temperature and
curvature (up the figure). Whereas the broadening with temperature is expected
since atomic displacement is a thermally-activated event, the broadening with
curvature is less so since the population and diffusivity of mobile surface
atoms are usually assumed to be curvature independent.

###  3.2 Diffusion and Activation Energy Results

The first part of Eq. 2.2 corresponds to the distribution of surface atom
displacements, with $D_{S}$ proportional to the variance and $n_{S}$ equal to
the probability mass of the distribution. Physically, $D_{S}$ is related to
the expected distance that a mobile surface atom travels in a given time,
whereas $n_{S}$ is related to the expected number of mobile surface atoms.

Figure 8: (a) Semi-log plot of surface diffusion coefficients $D_{S}$ for
several curvatures, found by fitting Eq. 2.2 to the displacements along
$\hat{\boldsymbol{z}}$. With $D_{S}$ as in Eq. 6, (b) the activation energy of
migration $\Delta H_{S}$ and (c) the surface diffusion coefficient prefactor
$D_{S}^{0}$ are shown as functions of curvature $\kappa_{1}$ for displacements
along $\hat{\boldsymbol{z}}$ and $\hat{\boldsymbol{\theta}}$. Error bars show
the standard deviation of the bootstrapping results and are negligible for
$D_{S}$.

Ref. | Method |  | Activation
---
Energy (eV)
| Prefactor
---
(Å2/ps)
[40] | Profile Decay | 0.91 |
[41] | Scratch method | 2.82 | $8.0\times 10^{10}$
[42] |  | Thin-film agglom-
---
erations
0.64 | $3.0\times 10^{-3}$
[43] | Scratch Method | 0.4 | $5.0\times 10^{-4}$
[44] | Crystalite profile | 0.35 | $5.0\times 10^{-5}$
[45] | Profile Decay | 1.1 | 5 to 10
[34] | MD simulation |  |
| (100) facet | 0.64, 0.84 | 8, 35
| (111) facet | 0.021, 0.038 | 2, 3.7
| (311) facet | 0.35, 0.42 | 8, 36
| (331) facet | 0.26, 0.34 | 6, 9.6
[33] | MD simulation |  |
| exchange | 0.15, 0.24 |
| jump | 0.11, 0.29, 0.34 |
| (100) facet | 0.41 |
[46] | MD simulation |  |
| (111) facet | 0.013, 0.014 | 0.45, 0.9

Table 1: A broad range of activation energies for surface diffusion have been
reported in the literature, though this could indicate that they are measuring
different populations of events [34]. The lower activation energies mostly
correspond to single adatom diffusion events.

The surface diffusion coefficients plotted in Fig. 8(a) were found by the
procedure described in Sec. 2.2, and depend on the assumption that the sets of
atoms belonging to the bulk and surface populations are fixed. Since $D_{S}$
is expected to follow the Arrhenius temperature dependence in Eq. 6, the
values are plotted with respect to $1/(k_{B}T)$ for each $\kappa_{1}$. The
activation energy for surface diffusion $\Delta H_{S}$ and the corresponding
prefactor $D_{S}^{0}$ were found from the slopes and the intercepts of the
plots of $D_{S}$ and are shown as functions of curvature in Figs. 8(b) and
(c). These reveal that $\Delta H_{S}$ appears to have a weak dependence on
$\kappa_{1}$, with the magnitudes of the slopes in Fig. 8(a) slightly
increasing with curvature. However, the apparent dependence on curvature is
likely a consequence of a finite expected lifetime of a mobile surface atom,
with the true surface diffusivity being independent of curvature as is
discussed in more detail below. The $\Delta H_{S}$ values found range from
$0.69$ to $0.79$ eV and are close to the median of the previously-reported
experimental values in Table 1. More specifically, the experimental values
range from $0.35$ to $2.82$ eV depending on the measurement method, and the
values from MD simulations range from $0.013$ to $0.84$ depending on the
specific migration event considered. It is significant that the activation
energy for any particular migration event in an MD simulation does not
adequately reflect the spectrum of events that occur in either an experimental
system or our MD simulations. Indeed, one benefit of the method developed in
Sec. 2.2 to measure surface diffusion is that the resulting diffusivities
naturally account for this spectrum of events.

The surface diffusion prefactor $D_{S}^{0}$ shown in Fig. 8(c) appears to be
an increasing function of curvature, but the dependence is small enough that
it would be difficult to measure experimentally. $D_{S}$ in Fig. 8(a) does
shift upward with increasing $\kappa_{1}$, though the trend is irregular with
the curves for $\kappa_{1}>0.0$ and $\kappa_{1}=0.0$ separated from those for
$\kappa_{1}<0.0$. While this suggests that $D_{S}^{0}$ (the vertical intercept
of the $D_{S}$ curves) could increase with curvature, the region plotted in
Fig. 8(a) is far enough from the vertical axis that the intercept value
depends more on the slope than the vertical offset. That is, the apparent
value of $D_{S}^{0}$ is highly sensitive to changes in $\Delta H_{S}$, and a
small systematic error in $\Delta H_{S}$ would be more than enough to account
for the apparent curvature dependence of $D_{S}^{0}$ in Fig. 8(c). The same
sensitivity occurs for experimental diffusion measurements as well, with the
surface diffusion prefactors in Table 1 increasing by 15 orders of magnitude
as the activation energies increase from $0.35$ to $2.82$ eV.

Figure 9: Semi-log plot of $n_{S}$ scaled by $V/(2Ar_{a})$ for several
curvatures (top), found by fitting Eq. 2.2 to the displacements along
$\hat{\boldsymbol{z}}$. With $n_{S}$ as in Eq. 7, the formation energy of a
mobile atom $E_{F}$ (bottom left) and the prefactor $\alpha_{0}$ (bottom
right) are shown as functions of curvature for displacements along
$\hat{\boldsymbol{z}}$ and $\hat{\boldsymbol{\theta}}$. The linear fits for
$E_{F}$ as a function of curvature are provided as guides for the eye. Error
bars show the standard deviation of the bootstrapping results and are
negligible for $n_{S}$. Figure 10: Cross-sectional view of a single Au
cylinder with a radius of $34$ Å simulated at $900$ K (left) and $1050$ K
(right). Atoms that have displaced more than $3$ Å after a $10$ ns simulation
are colored in yellow. For reference, $2.8$ Å is the approximate diameter of a
gold atom.

As described above, $n_{S}$ is the other significant quantity derived by
fitting atomic displacement data to Eq. 2.2, and is precisely defined as the
orientation-averaged fraction of atoms whose displacements along one of the
surface’s principal directions are characterized by the diffusivity $D_{S}$;
these atoms will be referred to as mobile surface atoms for brevity. This
fraction depends on the total number of atoms in the simulation though, which
for the punch hole geometry in particular is arbitrary. Hence, with $n_{S}$
modeled as in Eq. 7, $n_{S}V/(2r_{a}A)=\alpha_{0}\exp[-E_{F}/(k_{B}T)]$ is
plotted as a dimensionless quantity that roughly indicates the number of
effective mono-layers of mobile surface atoms. Note that the values for
$n_{S}V/(2r_{a}A)$ plotted in Fig. 9(a) are all of order one, consistent with
the expectation that roughly a single mono-layer of surface atoms are able to
diffuse by a short-circuit pathway. That said, the number of mobile surface
atoms can exceed the number of atoms in a single surface mono-layer at higher
temperatures; Figure 10 shows that one effective mono-layer of displaced atoms
at $900$ K increases to nearly two mono-layers at $1050$ K for this particular
crystallographic orientation. Notice that this does not necessarily indicate
that the concentration of mobile surface atoms is higher at $1050$ K than at
$900$ K at every instant in time, but rather that the population of atoms that
displaces on the surface at some point during the simulation (and therefore
contributes to $n_{S}$) increases with temperature. Indeed, the possibility of
a finite mobile surface atom lifetime is supported by the literature [32] and
by direct observations of atomic exchange between the two outermost effective
mono-layers in the simulation at $1050$ K shown in Fig. 10. This reinforces
that $n_{S}V/(2r_{a}A)$ should be considered as indicating the total number of
atoms that belonged to the mobile surface atom population at any point during
the simulation, rather than the surface concentration of atoms that are mobile
at every instant in time. The distinction is subtle but critical to an
accurate interpretation of our results since, for a fixed surface
displacement, the true diffusivity of a surface atom increases as the time
during which that atom belongs to the mobile surface population decreases.

The curves for $n_{S}V/(2r_{a}A)$ in Fig. 9(a) shift upward and decrease in
slope with increasing $\kappa_{1}$, appearing to all converge to a single
value that is independent of curvature at a temperature higher than $1050$ K
but much lower than the bulk melting point of gold at $1337$ K ($1440$ K for
the MEAM potential [38]). This would be quite surprising since the effect of
curvature on the mobile surface atom population would then be reversed at
still higher temperatures, with the apparent population of mobile surface
atoms being larger on negatively-curved surfaces than positively-curved ones
near the melting point. One possible resolution is that there could be surface
melting at a temperature below that of the bulk material [47], causing the
point of convergence to coincide with a change of transport mechanism and
making the apparent population of mobile surface atoms larger on positively-
curved solid surfaces at all temperatures. It is also significant that the
dependence of the slope on $\kappa_{1}$ makes the mobile surface atom
populations on positively-curved and negatively-curved surfaces increasingly
diverge at lower temperatures. Extrapolating the curves to $600$ K, the
temperature that is often used in thermal annealing experiments, suggests that
the apparent population of mobile surface atoms could be nearly twice as high
on a convex surface as on a concave one.

The slopes of the curves for $n_{S}V/(2r_{a}A)$ in Fig. 9(a) are precisely the
energy $E_{F}$ in Eq. 7, and are plotted with respect to $\kappa_{1}$ in Fig.
9(b). $E_{F}$ is interpreted as the characteristic energy required for an
immobile surface atom to become mobile. While we do not speculate as to the
mechanism, the magnitude of $E_{F}$ ($\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.27$ eV) is comparable to the single bond
energies of $0.15$ to $0.30$ eV routinely used in KMC simulations [30, 25].
Moreover, Fig. 9(b) shows that $E_{F}$ is a clearly decreasing function of
curvature given the magnitude of the random error. This is consistent with the
interpretation of $E_{F}$ above since the number of nearest neighbors of a
surface atom on average decreases with increasing curvature, reducing the
energetic barrier for an atom to join the population of mobile surface atoms
with increasing curvature. Notice that for the instantaneous concentration of
mobile surface atoms to be constant, the increased rate at which atoms join
the population of mobile surface atoms with increasing curvature would require
each atom to belong to that population for a shorter time, meaning that the
lifetime of a mobile surface atom would need to decrease with increasing
curvature.

A reasonable question is why identical reasoning should not apply to the
activation energy for surface diffusion $\Delta H_{S}$; if a decrease in
$E_{F}$ with curvature reflects a decrease in the average coordination number
of a surface atom, then one could expect $\Delta H_{S}$ to decrease with
curvature since the binding energy of the atoms to the surface would decrease
as well. One possible explanation relates to the proposed finite lifetime of a
mobile surface atom. If the rate of exchange of atoms between the bulk and
mobile surface atom populations increases with curvature, then the lifetime of
a mobile surface atom (i.e., the period during which the atom is able to
diffuse on the surface) would decrease with curvature. This could reduce the
apparent diffusivity that would be calculated given the assumption that the
set of atoms belonging to the mobile surface atom population is fixed, with
the result that the apparent activation energy of surface diffusion increases
as the mobile surface atom lifetime decreases. It is significant that this
could occur despite the average coordination number of a surface atom
decreasing with increasing curvature depending on the relative strengths of
the effects.

The quantity plotted in Fig. 9(c) is $\alpha_{0}$, corresponding to a
reference number of effective surface mono-layers in Eq. 7. The only physical
significance ascribed to this quantity is that
$\alpha_{0}\exp[-E_{F}/(k_{B}T)]$ should equal the number of effective mono-
layers of mobile surface atoms at the temperature at which the conjectured
surface melting occurs. While $\alpha_{0}$ does appear to be a decreasing
function of curvature, the dependence is not believed to be significant, with
a small systematic error in $E_{F}$ being more than enough to account for the
curvature dependence of $\alpha_{0}$ in Fig. 9(c).

It is worthwhile commenting briefly on the two values plotted for the $z$
direction in Figs. 8(b), 8(c), 9(b) and 9(c) for $\kappa_{1}=0$. There are no
results for the $\theta$ direction since this direction is poorly defined in
the slab geometry. Instead, the $\hat{\boldsymbol{x}}$ displacements were
aggregated into a separate data set, with the derived values colored the same
way as for $\hat{\boldsymbol{z}}$ displacements since the results along the
$z$ and $x$ directions should coincide given sufficient sampling of crystal
orientations. As for the simulations of cylinders and punches, the difference
between diffusion along the $\theta$ and $z$ directions (corresponding to
$\hat{\kappa}_{1}$ and $\hat{\kappa}_{2}$ respectively) is subtle. It appears
in Figure 8(b) that $\Delta H_{S}$ is slightly higher for diffusion in the $z$
direction; this is attributed to the relatively small number of crystal
orientations sampled, since the magnitude of the difference is comparable to
that between the $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{z}}$
displacements for the slab geometry.

Transport theory indicates that the rate of mass transport by surface
diffusion should be directly proportional to the product of the number of
mobile surface atoms and the surface diffusivity. Given the Arrhenius
dependencies of Eqs. 6 and 7, this implies that the rate of mass transport by
surface diffusion should have an Arrhenius dependence with an activation
energy of $E_{F}+\Delta H_{S}$ and a prefactor proportional to
$\alpha_{0}D^{0}_{S}$. These quantities are plotted in Fig. 11 as functions of
the curvature, and appear to be independent of curvature within the estimated
uncertainty. While it is initially surprising that the curvature dependencies
of $D_{S}$ and $n_{S}$ should be nearly equal and opposite, decades of
literature has assumed that mass transport by surface diffusion is independent
of curvature [29]. This means that the real question is why the analysis
performed here found $D_{S}$ and $n_{S}$ to have any curvature dependence at
all.

Figure 11:  The rate of mass transport by surface diffusion is expected to
have an Arrhenius dependence with an activation energy $E_{F}+\Delta H_{S}$
and a prefactor proportional to $\alpha_{0}D^{0}_{S}$. These quantities are
independent of curvature to within the uncertainty estimates, despite the
apparent curvature dependencies in Figs. 8 and 9.

As suggested in Sec. 2.2, we propose that the apparent dependence on curvature
is not introduced by the simulation setup, but by one of the assumptions
behind Eq. 2.2. Specifically, Eq. 2.2 assumes that the set of atoms belonging
to the mobile surface atom population is fixed for the entire simulation
duration, and that no atoms are exchanged with the immobile bulk atom
population. If this were true, then the variance of the displacement
distribution of the fixed population of mobile surface atoms would indeed
increase as $tD_{S}$ as is assumed above. However, if only the concentration
of mobile surface atoms is fixed, then individual atoms could be exchanged
between the mobile surface atom and immobile bulk populations at a curvature-
dependent rate over the simulation duration. Increasing this rate would
decrease the mobile surface atom lifetime and, as the simulation duration
increases relative to this lifetime, have two significant effects on derived
quantities. First, the apparent population of mobile surface atoms would
increase since many atoms would join and leave the population of mobile
surface atoms at some point during the simulation. Second, the apparent
surface diffusivity would decrease since any single mobile surface atom would
displace on the surface for a decreasing fraction of the simulation. Figure 12
shows that this appears to be the case, with $n_{S}V/(2r_{a}A)$ on the left
increasing and $D_{S}$ in the middle decreasing with the simulation duration
$\Delta t$. The reasoning in Sec. 2.2 further suggests that the flux of
surface atoms and the product $n_{S}D_{S}$ should be independent of $\Delta
t$, and this is also supported by the plot on the right of Fig. 12 after an
initial transient related to the characteristic time for an atomic jump.
Notice that $n_{S}$ and $D_{S}$ are expected to approach asymptotic values for
simulation durations much longer than the mobile surface atom lifetime when
individual atoms can join and leave the mobile surface atom population
multiple times, but a direct investigation is not computationally feasible at
present.

Figure 12: The values plotted are calculated from a data set made up of atomic
displacements found after time $\Delta t$. The data set sampled from for each
$\Delta t$ is made up from four different simulations of a 30 Å cylinder at
$1050$ $K$ each with different crystal orientations. Mobile atom exchange with
an immobile atom population increases with time and this is reflected in the
calculated values for $D_{S}$ and $n_{S}$ though their product is relatively
constant.

The apparent curvature dependence of $D_{S}$ and $n_{S}$ can therefore be
explained as a consequence of an increasing rate of exchange between the
mobile surface atom and immobile bulk atom populations with increasing
curvature. Notice that the overall rate of mass transport by surface diffusion
does not need to depend on this rate of exchange; the instantaneous
concentration of mobile surface atoms and the rate at which they jump to
neighboring surface sites could easily be independent of both the rate of
exchange and the surface curvature. Unfortunately, there does not seem to be
any way to directly measure these quantities by means of MD simulations, at
least not without explicitly defining what is meant by a mobile surface atom,
and (as already described in Sec. 2.2) it is not obvious whether this can be
done in an unbiased way. Our conclusion is that this unexpected phenomenon
could make the apparent surface diffusivity reported by simulation studies
dependant on not only the surface curvature but on the simulation time
(relative to the mobile surface atom lifetime). Moreover, making the standard
assumption that the apparent population of mobile surface atoms is constant
and independent of curvature and temperature could result in the activation
energy for mass transport by surface diffusion being reported as $\Delta
H_{S}$ alone, differing by $E_{F}$ ($\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.26$ eV) from the value reported here.

##  4 Conclusion

Many of the desirable properties of nanoporous metals in general, and np-Au in
particular, are a consequence of an exceptionally high specific surface area
$S_{A}$. This high $S_{A}$ provides the driving force for the coarsening of
the microstructure at moderate temperatures by surface diffusion, though
directly measuring mass transport rates along individual ligaments by in-situ
imaging is challenging due to the small features and short times involved. An
alternative approach is to infer the surface diffusivity by means of a scaling
argument, assuming that the microstructure coarsens in a self-similar way.
However, there is no consensus in the literature that the microstructure
coarsening is actually self-similar, and the activation energies found by this
procedure can be as low as $0.353$ eV [31] or as high as $1.33$ eV [1]. This
results in an uncertainty in the rate of mass transport by surface diffusion
that makes it difficult to precisely predict the microstructural change that
would result from a given thermal annealing schedule, or more generally, to
construct processing routes to a nanoporous material with a desired
microstructure.

Given this situation, atomistic simulations are a natural candidate to measure
the rates of mass transport by surface diffusion in a simulated material.
Connecting such simulation results to the quantities that appear in classical
transport theory is not trivial though, with the definition of various
quantities like the chemical potential remaining unclear at the atomic scale.
The approach developed here compares the distribution of atomic displacements
in the absence of a chemical potential gradient to that expected for an
anisotropic random walker. The resulting values for the surface diffusivity
and mobile surface atom population appear to suggest subtle dependencies on
the surface mean curvature, despite the overall rate of mass transport by
surface diffusion being curvature independent. It is proposed that the
apparent curvature dependence of the surface diffusivity and mobile surface
atom population is a consequence of the surface atoms having a finite lifetime
during which they remain mobile (a possibility supported by the literature
[32]), after which they are exchanged with one of the population of
effectively immobile bulk atoms. This could occur in such a way that the
concentration of mobile surface atoms remains constant, making the overall
rate of mass transport by surface diffusion independent of the rate of
exchange. While an increasing rate of exchange with increasing surface
curvatures would be able to explain our results, directly quantifying the rate
of exchange would require a precise definition of what constitutes a mobile
surface atom, and is not performed here.

The rate of mass transport is directly proportional to the product of the
surface diffusivity and the population of mobile surface atoms, and is found
in Fig. 11 to have an activation energy of $\raise
0.73193pt\hbox{$\scriptstyle\sim$}1.01$ eV, close to the median of the
experimentally-determined activation energies reported in Table 1. This
includes an activation energy of $\raise
0.73193pt\hbox{$\scriptstyle\sim$}0.27$ eV related to the activation energy
for the exchange mechanism and the apparent dependence of the mobile surface
atom population on temperature. Given that the concentration of mobile surface
atoms is conjectured to be constant, it would be easy for a researcher
studying surface diffusion with MD simulations to assume that the population
of mobile surface atoms should also be constant, and that the rate of mass
transport should depend on temperature only through the surface diffusivity.
Measuring the surface diffusivity by comparing the displacement distribution
of surface atoms to that of an anisotropic random walker would then lead to
the the activation energy for mass transport by surface diffusion being
incorrectly reported as only $\raise 0.73193pt\hbox{$\scriptstyle\sim$}0.74$
eV, a difference that would significantly change the predicted thermal
annealing schedules required to achieve a particular ligament dimension in
experiments.

Finally, it is significant that the activation energy for mass transport found
here is appreciably higher than that for surface diffusion for any of the MD
simulations reported in Table 1. Overall, this study underlines the importance
of considering the possible effects of a range of migration events and the
finite lifetime of a mobile surface atom when extracting surface diffusion
values in simulated materials.

## Data Availability

The raw and processed data required to reproduce these findings are available
for download in the supplementary materials.

## Acknowledgements

C.M.E, S.M.S., E.S and J.K.M. were supported by the National Science
Foundation under Grant No. 2003849.

## References

  * [1] I. McCue, J. Stuckner, M. Murayama, M. J. Demkowicz, Gaining new insights into nanoporous gold by mining and analysis of published images, Scientific Reports 8 (1) (Apr. 2018).  doi:10.1038/s41598-018-25122-3.
  * [2] A. Wittstock, J. Biener, M. Bäumer, Nanoporous gold: a new material for catalytic and sensor applications, Physical Chemistry Chemical Physics 12 (40) (2010) 12919.  doi:10.1039/c0cp00757a.
  * [3] J. Biener, M. M. Biener, R. J. Madix, C. M. Friend, Nanoporous gold: Understanding the origin of the reactivity of a 21st century catalyst made by pre-columbian technology, ACS Catalysis 5 (11) (2015) 6263–6270.  doi:10.1021/acscatal.5b01586.
  * [4] D. Wang, Z. Hu, G. Peng, Y. Yin, Surface energy of curved surface based on lennard-jones potential, Nanomaterials 11 (3) (2021) 686.  doi:10.3390/nano11030686.
  * [5] J. Qi, J. Zeng, F. Zhao, S. H. Lin, B. Raja, U. Strych, R. C. Willson, W.-C. Shih, Label-free, in situ SERS monitoring of individual DNA hybridization in microfluidics, Nanoscale 6 (15) (2014) 8521–8526.  doi:10.1039/c4nr01951b.
  * [6] L. Zhang, H. Chang, A. Hirata, H. Wu, Q.-K. Xue, M. Chen, Nanoporous gold based optical sensor for sub-ppt detection of mercury ions, ACS Nano 7 (5) (2013) 4595–4600.  doi:10.1021/nn4013737.
  * [7] F. Ruffino, M. G. Grimaldi, Nanoporous gold-based sensing, Coatings 10 (9) (2020) 899.  doi:10.3390/coatings10090899.
  * [8] C. A. R. Chapman, H. Chen, M. Stamou, J. Biener, M. M. Biener, P. J. Lein, E. Seker, Nanoporous gold as a neural interface coating: Effects of topography, surface chemistry, and feature size, ACS Applied Materials & Interfaces 7 (13) (2015) 7093–7100.  doi:10.1021/acsami.5b00410.
  * [9] P. Daggumati, S. Appelt, Z. Matharu, M. L. Marco, E. Seker, Sequence-specific electrical purification of nucleic acids with nanoporous gold electrodes, Journal of the American Chemical Society 138 (24) (2016) 7711–7717.  doi:10.1021/jacs.6b03563.
  * [10] E. Şeker, W.-C. Shih, K. J. Stine, Nanoporous metals by alloy corrosion: Bioanalytical and biomedical applications, MRS bulletin 43 (1) (2018) 49–56.
  * [11] Y. Yu, L. Gu, X. Lang, C. Zhu, T. Fujita, M. Chen, J. Maier, Li storage in 3d nanoporous au-supported nanocrystalline tin, Advanced Materials 23 (21) (2011) 2443–2447.  doi:10.1002/adma.201004331.
  * [12] J. M. Gonçalves, A. Kumar, M. I. da Silva, H. E. Toma, P. R. Martins, K. Araki, M. Bertotti, L. Angnes, Nanoporous gold-based materials for electrochemical energy storage and conversion, Energy Technology 9 (5) (2021) 2000927\.  doi:10.1002/ente.202000927.
  * [13] J. Erlebacher, M. J. Aziz, A. Karma, N. Dimitrov, K. Sieradzki, Evolution of nanoporosity in dealloying, Nature 410 (6827) (2001) 450–453.  doi:10.1038/35068529.
  * [14] J. Snyder, P. Asanithi, A. B. Dalton, J. Erlebacher, Stabilized nanoporous metals by dealloying ternary alloy precursors, Advanced Materials 20 (24) (2008) 4883–4886.  doi:10.1002/adma.200702760.
  * [15] J. Weissmüller, R. C. Newman, H.-J. Jin, A. M. Hodge, J. W. Kysar, Nanoporous metals by alloy corrosion: Formation and mechanical properties, MRS Bulletin 34 (8) (2009) 577–586.  doi:10.1557/mrs2009.157.
  * [16] J. Erlebacher, R. Seshadri, Hard materials with tunable porosity, MRS Bulletin 34 (8) (2009) 561–568.  doi:10.1557/mrs2009.155.
  * [17] T. Fujita, L.-H. Qian, K. Inoke, J. Erlebacher, M.-W. Chen, Three-dimensional morphology of nanoporous gold, Applied Physics Letters 92 (25) (2008) 251902.  doi:10.1063/1.2948902.
  * [18] M. Ziehmer, K. Hu, K. Wang, E. T. Lilleodden, A principle curvatures analysis of the isothermal evolution of nanoporous gold: Quantifying the characteristic length-scales, Acta Materialia 120 (2016) 24–31.  doi:10.1016/j.actamat.2016.08.028.
  * [19] E. Şeker, J. T. Gaskins, H. Bart-Smith, J. Zhu, M. L. Reed, G. Zangari, R. Kelly, M. R. Begley, The effects of post-fabrication annealing on the mechanical properties of freestanding nanoporous gold structures, Acta Materialia 55 (14) (2007) 4593–4602.  doi:10.1016/j.actamat.2007.03.018.
  * [20] D. M. Turner, S. R. Niezgoda, S. R. Kalidindi, Efficient computation of the angularly resolved chord length distributions and lineal path functions in large microstructure datasets, Modelling and Simulation in Materials Science and Engineering 24 (7) (2016) 075002.  doi:10.1088/0965-0393/24/7/075002.
  * [21] P. Gao, Z. Lv, Y. Song, M. Song, P. Qian, Evolution of morphology and microstructure of coarsened nanoporous gold studied by automatic thresholding and image recognition algorithms, Scripta Materialia 226 (2023) 115256.  doi:10.1016/j.scriptamat.2022.115256.
  * [22] Y. chen Karen Chen, Y. S. Chu, J. Yi, I. McNulty, Q. Shen, P. W. Voorhees, D. C. Dunand, Morphological and topological analysis of coarsened nanoporous gold by x-ray nanotomography, Applied Physics Letters 96 (4) (2010) 043122.  doi:10.1063/1.3285175.
  * [23] Y. K. Chen-Wiegart, S. Wang, Y. S. Chu, W. Liu, I. McNulty, P. W. Voorhees, D. C. Dunand, Structural evolution of nanoporous gold during thermal coarsening, Acta Materialia 60 (12) (2012) 4972–4981.  doi:10.1016/j.actamat.2012.05.012.
  * [24] J. Erlebacher, An atomistic description of dealloying, Journal of The Electrochemical Society 151 (10) (2004) C614.  doi:10.1149/1.1784820.
  * [25] Y. Li, B.-N. D. Ngô, J. Markmann, J. Weissmüller, Topology evolution during coarsening of nanoscale metal network structures, Physical Review Materials 3 (7) (Jul. 2019).  doi:10.1103/physrevmaterials.3.076001.
  * [26] K. Kolluri, M. J. Demkowicz, Coarsening by network restructuring in model nanoporous gold, Acta Materialia 59 (20) (2011) 7645–7653.  doi:10.1016/j.actamat.2011.08.037.
  * [27] M. Hakamada, M. Mabuchi, Thermal coarsening of nanoporous gold: Melting or recrystallization, Journal of Materials Research 24 (2) (2009) 301–304.  doi:10.1557/jmr.2009.0037.
  * [28] A. Chen, S. Shi, F. Liu, Y. Wang, X. Li, J. Gu, X. Xie, Effect of annealing atmosphere on the thermal coarsening of nanoporous gold films, Applied Surface Science 355 (2015) 133–138.  doi:10.1016/j.apsusc.2015.07.065.
  * [29] C. Herring, Effect of change of scale on sintering phenomena, Journal of Applied Physics 21 (4) (1950) 301–303.  doi:10.1063/1.1699658.
  * [30] J. Erlebacher, Mechanism of coarsening and bubble formation in high-genus nanoporous metals, Physical Review Letters 106 (22) (Jun. 2011).  doi:10.1103/physrevlett.106.225504.
  * [31] H. Jeon, N.-R. Kang, E.-J. Gwak, J. il Jang, H. N. Han, J. Y. Hwang, S. Lee, J.-Y. Kim, Self-similarity in the structure of coarsened nanoporous gold, Scripta Materialia 137 (2017) 46–49.  doi:10.1016/j.scriptamat.2017.05.009.
  * [32] G. Antczak, G. Ehrlich, Jump processes in surface diffusion, Surface science reports 62 (2) (2007) 39–61.
  * [33] F. Baletto, C. Mottet, R. Ferrando, Molecular dynamics simulations of surface diffusion and growth on silver and gold clusters, Surface Science 446 (1-2) (2000) 31–45.  doi:10.1016/s0039-6028(99)01058-4.
  * [34] C. Liu, J. Cohen, J. Adams, A. Voter, EAM study of surface self-diffusion of single adatoms of fcc metals ni, cu, al, ag, au, pd, and pt, Surface Science 253 (1-3) (1991) 334–344.  doi:10.1016/0039-6028(91)90604-q.
  * [35] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. t. in Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271 (2022) 108171.  doi:10.1016/j.cpc.2021.108171.
  * [36] G. Grochola, S. P. Russo, I. K. Snook, On fitting a gold embedded atom method potential using the force matching method, The Journal of Chemical Physics 123 (20) (2005) 204719.  doi:10.1063/1.2124667.
  * [37] C. A. Becker, Z. T. Trautt, L. M. Hale, 2005–grochola-g-russo-s-p-snook-i-k–au (Feb 2011).
  * [38] B.-J. Lee, J.-H. Shim, M. I. Baskes, Semiempirical atomic potentials for the fcc metals cu, ag, au, ni, pd, pt, al, and pb based on first and second nearest-neighbor modified embedded atom method, Physical Review B 68 (14) (Oct. 2003).  doi:10.1103/physrevb.68.144112.
  * [39] R. P. Brent, Algorithms for minimization without derivatives, Courier Corporation, 2013.
  * [40] T.-S. Lin, Y.-W. Chung, Measurement of the activation energy for surface diffusion in gold by scanning tunneling microscopy, Surface Science 207 (2-3) (1989) 539–546.  doi:10.1016/0039-6028(89)90140-4.
  * [41] M. McLean, J. Hirth, Surface self-diffusion on gold, Surface Science 12 (2) (1968) 177–188.  doi:10.1016/0039-6028(68)90122-2.
  * [42] I. Beszeda, I. Szabó, E. Gontier-Moya, Morphological evolution of thin gold films studied by auger electron spectroscopy in beading conditions, Applied Physics A 78 (7) (2004) 1079–1084.  doi:10.1007/s00339-003-2170-0.
  * [43] H. Göbel, P. von Blanckenhagen, A study of surface diffusion on gold with an atomic force microscope, Surface Science 331-333 (1995) 885–890.  doi:10.1016/0039-6028(95)00206-5.
  * [44] M. Drechsler, J. Métois, J. Heyraud, Surface self-diffusion studied by microscopic measurements of crystallite profile evolutions, Surface Science 108 (3) (1981) 549–560.  doi:10.1016/0039-6028(81)90564-1.
  * [45] D. Porath, E. Bar-Sadeh, M. Wolovelsky, A. Grayevsky, Y. Goldstein, O. Millo, Annealing study of gold films using scanning tunneling microscopy, Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 13 (3) (1995) 1165–1170.  doi:10.1116/1.579855.
  * [46] G. Boisvert, L. J. Lewis, Single-atom diffusion on silver and gold surfaces, MRS Proceedings 317 (1993).  doi:10.1557/proc-317-71.
  * [47] G. Guenther, O. Guillon, Models of size-dependent nanoparticle melting tested on gold, Journal of materials science 49 (2014) 7915–7932.